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Abstract 

Micro-air  vehicles  (MAV)s  provide  a  valuable  and  low  observable  way  to  do  the 
jobs  that  the  Air  Force  deems  to  be  dull,  dirty  and  dangerous.  Basing  the  design  of  an 
MAV  wing  on  that  of  a  biological  counterpart  will  provide  a  proven  design  that  is 
capable  of  achieving  the  mission  requirements.  This  research  is  designed  to  analyze  the 
design  and  manufacturing  of  a  wing  based  off  the  Manduca  Sexta. 

Inaccuracies  in  the  fiber  orientation  can  result  in  substantial  changes  in  the 
material  properties.  Experimental  vibration  data  of  composite  material  samples 
manufactured  using  a  three-ply  [0/90/0]  small  non-homogenous  fiber  composite  provided 
results  that  varied  over  33  percent  from  analytical  results.  Since  the  material  was  to  be 
used  in  the  manufacturing  of  a  biologically  inspired  MAV,  it  was  important  to  understand 
the  cause  of  the  variance  in  the  measured  material  properties  so  that  they  could  be  taken 
into  account  for  the  design  and  manufacturing  of  the  MAV  wing. 

An  analysis  was  performed  on  the  material  to  verify  that  it  matched  specified 
material  properties.  Inaccuracies  in  the  manufacturing  of  the  composite  samples  were 
taken  into  account;  specifically  ply  orientation,  cut  angle,  and  material  thickness  were 
examined.  Using  finite  element  analysis  (FEA),  it  was  determined  that  a  misalignment  in 
fiber  orientation  of  less  than  five  degrees  combined  with  resulting  short  fiber  effects 
accounts  for  the  difference  between  analytical  and  experimental  results.  Using  an  optical 
microscope,  variances  in  the  ply  orientation  was  observed  confirming  the  FEA  results. 
Possible  inaccuracies  in  the  composite  material  were  taken  into  consideration  during  the 
design  and  construction  of  the  MAV  wing. 


v 


A  FEA  model  of  the  engineered  MAV  wing  was  developed  with  the  carbon-fiber 
composites  inaccuracies  in  mind.  To  allow  for  changes  to  the  model  to  be  made  quickly, 
the  FEA  model  was  generated  using  a  developed  MATLAB  code  that  generated  finite 
element  input  files  to  be  solved  using  ABAQUS,  a  finite  element  program.  The 
developed  MATLAB  code  generated  beam  cross-sections  for  the  composite  material 
elements  based  upon  the  input  ply  orientations  and  inaccuracies  and  assigned  corrected 
densities  of  each  of  the  beam  elements.  It  also  allowed  for  idiosyncrasies  of  the 
composite  vein  structure  of  the  wing  to  easily  be  changed  and  evaluated.  Since  multiple 
FEA  input  files  could  be  generated  quickly,  an  analysis  was  performed  in  order  to 
detennine  the  effects  of  the  misalignment  of  the  ply  orientation  angles  on  the  MAV  wing 
model. 
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A  FINITE  ELEMENT  ANALYSIS  OF  A  CARBON  FIBER  COMPOSITE  MICRO 

AIR  VEHICLE  WING 

I.  Introduction 


1.1.  Research  Objective 

Current  research  done  by  O’Hara  [1]  on  the  Manduca  Sexta  has  provided  a 
preliminary  design  for  a  manufactured  bio-inspired  Micro  Air  Vehicle  (MAV)  wing 
based  upon  the  structural  characteristics  of  the  forewing  of  the  Manduca  Sexta.  This 
preliminary  design  consists  of  a  unidirectional  carbon  fiber  composite  vein  structure  with 
a  Kapton  membrane  manufactured  at  the  Air  Force  Institute  of  Technology  (AFIT),  see 
Figure  1.1.  The  materials  for  the  wing  were  chosen  based  upon  the  need  of  the  MAV 
wing  to  have  both  high  strength  and  low  density,  successful  testing  compared  to  metallic 
materials  and  previous  research  found  showing  success  of  carbon  fiber  wings  in  other 
MAV  applications  [1]  [2]  [3].  The  veins  of  the  wing  were  cut  out  of  a  single  sheet  of 
[0/90/0]  three-ply  unidirectional  carbon-fiber  composite,  based  on  the  geometry  of  the 
Manduca  Sexta,  and  a  “best  guess”  method  of  determine  the  vein  width  in  order  to  create 
a  wing  that  could  withstand  the  stresses  and  strains  of  testing. 


Figure  1.1:  Manufactured  MAV  Wing  Developed  by  O'Hara  [2], 
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While  work  done  by  DeLeon  has  shown  that  the  design  of  O’Hara’s  wing  shares 
some  characteristics  of  the  hawkmoth  wing,  some  further  understanding  is  necessary  to 
create  a  manufactured  wing  that  will  more  closely  match  the  structural  characteristics  of 
the  Manduca  Sexta  [1].  It  is  important  that  any  design  changes  be  made  with  a  full 
understanding  of  the  structure  at  hand,  so  that  a  better  design  will  result.  This  will  not 
only  lead  to  a  better-designed  wing,  but  also  a  more  efficient  design  process. 

The  carbon  fiber  composite  used  in  the  engineered  wing  allowed  for  a  stiffer, 
lightweight  design.  However,  there  were  concerns  over  the  effectiveness  of  the  carbon 
fiber  composite  in  the  wing  design.  Would  the  carbon  fiber  composite  retain  its  material 
properties  given  the  small  size  of  the  veins?  How  does  the  fact  that  the  veins  are  laser  cut 
affect  the  structure  of  the  engineered  wing?  How  does  the  anisotropic  nature  of  the 
composite  affect  the  wing? 

These  questions  are  worth  looking  at  if  the  use  of  the  three-ply  unidirectional 
composite  in  the  MAV  wing  is  to  be  considered.  Therefore,  the  objective  of  this  thesis  is 
to  perfonn  an  analysis  of  the  engineered  MAV  wing  structure  and  the  carbon  fiber 
composite  material  of  which  it  is  mainly  composed.  Furthermore,  the  generated  finite 
element  model  should  be  easily  modified  so  that  the  effects  of  design  or  material  changes 
to  the  wing  can  easily  be  made.  This  will  allow  for  future  designs  to  be  analyzed  prior  to 
construction,  thus  saving  considerable  time  in  the  design  process.  In  performing  the 
structural  analysis  using  the  finite  element  method,  many  of  the  questions  above  can  be 
answered,  and  a  better  understanding  of  the  manufactured  MAV  wing  can  be 
accomplished. 
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1.2.  Background 

The  demand  for  unmanned  intelligence,  surveillance  and  reconnaissance  (ISR)  has  grown 


significantly  in  recent  decades.  Unmanned  aerial  vehicles  (UAVs)  are  unrestricted  by  a 
pilot’s  physical  limitations  and  can  provide  unprecedented  loiter  times,  range,  and  cost 
effectiveness.  MAVs  have  been  proposed  by  many  due  to  low  cost,  tremendous 
maneuverability  and  inconspicuous  operation,  and  are  seen  by  many  as  the  next 
revolution  in  the  field  of  UAVs.  The  Defense  Advanced  Research  Agency’s  (DARPA’s) 
current  vision  for  the  optimal  MAY  is  outlined  in  Table  1.1. 


Table  1.1:  MAY  Design  Requirements  [4]. 


Specification 

Requirements 

Details 

Size 

<15.24  cm 

Maximum  Dimension 

Weight 

~100g 

Objective  GTOW 

Range 

1  to  10  km 

Operational  Range 

Endurance 

60  min 

Loiter  time  on  station 

Altitude 

<150  m 

Operational  ceiling 

Speed 

15  m/s 

Maximum  flight  speed 

Payload 

20  g 

Mission  dependent 

Cost 

$1,500 

Maximum  cost,  2009  USD 

A  quick  glance  at  the  requirements  will  reveal  a  significant  difference  between 
current  aircraft  requirements  and  the  requirements  of  an  MAV.  With  such  contrast  in 
vehicle  requirements,  it  is  only  natural  that  the  design  of  a  vehicle  will  be  different. 
Using  nature  as  a  guide,  there  are  many  examples  of  insects  that  use  flapping  wings  as 
highly  efficient  mechanisms  for  flight  [5].  The  use  of  such  a  flapping  mechanism,  where 
the  wing  would  produce  thrust  as  well  as  generate  lift,  represents  a  dramatic  shift  from 
the  traditional  quasi-rigid  wing  or  rotorcraft  with  a  separate  propulsion  system.  Research 
of  the  countless  number  of  flying  insects  present  in  the  world  provide  examples  of 
structures  that  are  structurally  capable  of  flight  and  have  the  potential  to  meet  DARPA’s 
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requirements  for  MAVs.  These  structures  would  provide  a  valid  aerodynamic  model  for 
research.  This  would  reduce  the  significant  amount  of  time  associated  with  solving  the 
low  Reynolds’s  number  and  unsteady  aerodynamics  of  a  flapping  MAV.  One  insect  that 
possesses  the  capabilities  to  meet  these  requirements  is  the  Manduca  Sexta,  also  known 
as  the  North  American  Hawkmoth  or  simply  as  the  hawkmoth  (see  Figure  1.2)  [6]. 
Insects  such  at  the  Manduca  Sexta  have  the  benefit  of  having  undergone  millions  of  years 
of  evolution  and  have  evolved  into  highly  specialized  and  efficient  systems.  With  such  a 
valid  aerodynamic  model,  if  the  structural  properties  of  the  hawkmoth  could  be  matched, 
a  capable  MAV  could  be  developed. 


Figure  1.2:  An  Adult,  Female  Manduca  Sexta  (Flawkmoth). 


The  Manduca  Sexta  forewing  was  selected  as  the  basis  for  the  bio-inspired 
design.  The  forewing,  which  produces  most  of  the  lift  for  the  moth,  has  been  selected 
and  studied  by  Norris,  and  has  been  determined  to  be  an  ideal  candidate  for  a  wing  [6]. 
O’Hara  has  continued  the  research  of  Norris  and  has  established  general  material 
properties  and  characteristic  dimensions  for  the  hawkmoth  wing,  and  as  developed  a 
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preliminary  design  for  an  MAV  wing  based  on  the  his  findings,  using  carbon-fiber  and 
Kapton  (see  Figure  1.1)  [2],  O’Hara’s  research  provides  the  basis  for  a  bio-inspired 
MAV  wing,  and  for  the  objective  of  this  thesis  [2]. 

1.3.  Motivation 

As  war  fighters  are  taxed  with  ever  increasing  difficult  situations,  it  is  crucial  that 
they  be  made  aware  of  what  is  to  come;  whether  it  be  over  a  hill,  in  a  cave,  or  in  a  room 
of  an  unknown  building  [7],  DARPA’s  vision  for  a  small,  lightweight  MAV  will  allow 
war  fighters  to  observe  the  situation  while  the  vehicles  small  size  will  allow  it  to  blend  in 
and  be  nearly  indistinguishable  from  the  natural  insect  population.  The  capability  for  a 
vehicle  to  be  able  to  observe  a  situation  while  hiding  in  plain  sight  allows  it  the  ability  to 
gather  ISR  previously  unobtainable  by  traditional  UAVs. 

The  ability  for  a  vehicle  to  successfully  navigate  buildings  requires  it  to  be  highly 
maneuverable.  This  is  quite  difficult  to  achieve  with  traditional  quasi-rigid  winged 
vehicles,  since  forward  motion  is  constantly  required  to  maintain  lift.  It  can  be  said  that 
in  order  to  meet  such  a  challenging  design  criteria  it  is  best  to  turn  to  nature.  Norris  et  al. 
stated  that  current  research  finds  a  vested  interest  in  this  flapping  wing  MAV  (FWMAV), 
and  that  science  is  compelled  to  mimic  the  elegant  (and  efficient)  designs  that  nature  has 
developed  for  its  flapping  wing  design  [8]. 

Current  designs  of  MAVs,  shown  in  Figure  1.3,  show  that  much  research  is  still 
required  in  order  to  achieve  a  design  that  is  capable  to  fulfill  DARPA’s  design  criteria. 
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Micro  Commercial  Allied  Aerospace  iSTAR 

Rc  Heli  [350g  /  15min]  [1500g/15  min  9  in  dia] 


Mentor  [550g/10  min] 


MicroSTAR  [110g/25min] 
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Objective 

[100g/60min] 


Microbat  [I0g/15  min] 


Endurance  (min) 


Black  WidOW [80g/22min] 


Figure  1.3:  View  of  Recent  MAY  Developments  [9]. 


1.4.  The  Concept  of  Biological  Inspiration 

Looking  to  nature’s  countless  examples  of  flying  creatures  is  not  a  new  concept. 

The  Wright  Brothers  built  and  tested  airfoil  designs  based  upon  the  shapes  of  the  wings 
of  birds,  eventually  leading  to  the  creation  of  the  Wright  Flyer.  It  is  easy,  however,  for 
scientists  to  forget  about  nature’s  influence  on  flight  and  focus  energy  into  determining 
the  ‘best  way’  for  something  to  work,  as  opposed  to  examining  what  is  ‘known’  to  work 
in  nature  and  determining  how  it  works.  Without  looking  at  everything  that  is  currently 
known,  many  find  themselves  effectively  reinventing  the  wheel  instead  of  using 
knowledge  and  information  already  known. 
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Biologically  inspired  flight  can  be  said  to  be  the  foundation  of  flapping  wing 
flight.  Erich  von  Holst  is  considered  an  early  pioneer  in  the  study  of  bio-inspired 
flapping  wing  vehicles.  He  and  his  colleagues  studied  the  biological  and  aerodynamics 
of  flying  animals  in  the  1930’s  and  1940’s,  even  developing  several  successful  flapping 
vehicles  based  on  the  study  of  birds  and  insects  (See  Figure  1.4.  [10]). 


Figure  1.4:  Erich  von  Holst  with  a  Flapping  Air  Vehicle  Based  upon  a  Swan  [10]. 


In  the  period  following  World  War  II,  it  seems  that  the  importance  of  bio-inspired 
work  was  forgotten  until  only  recently.  Aaron  Norris  found  that  there  was  a  significant 
lack  of  literature  present  on  the  area  of  bio-inspired  flight  [6].  This  especially  held  true  in 
the  area  of  insects  that  would  likely  be  mimicked  in  order  to  create  an  MAY  based  upon 
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DARPA’s  criteria.  Furthermore,  Norris  and  DeLeon  both  found  that  the  majority  of 
literature  that  exists  dealt  with  the  aerodynamics  associated  with  MAV,  and  that  there 
was  a  significant  lack  of  information  about  the  structures  of  MAVs  [6]  (See  Figure  1.5. 
[1]).  In  order  to  create  a  feasible  MAV  wing,  more  research  was  needed. 


Figure  1.5:  Current  Research  Intrest  of  Flapping  Flight  [6]. 

1.5.  A  Look  at  the  Manduca  Sexta 

It  is  important  to  look  at  some  of  the  previous  work  performed  on  the  Manduca 
Sexta  to  better  understand  what  is  trying  to  be  matched  for  the  creation  of  a  flapping 
micro  air  vehicle  (FMAV).  Looking  at  previous  work  done  [1]  [2],  one  will  better 
understand  the  motivation  behind  the  project,  as  well  as  the  objectives  for  an  engineered 


wing.  This  section  will  primarily  focus  on  the  work  done  by  Norris  and  O’Hara,  its 
relation  to  the  Manduca  Sexta  and  application  for  FWMAVs. 

Norris  studied  the  forewing  of  the  Manduca  Sexta,  selecting  the  hawkmoth 
because  it  was  readily  available  for  study  and  could  meet  the  requirements  set  forth  by 
DARPA  for  an  MAV.  Norris  “liberated”  (i.e.  separated)  the  wings  from  a  large  sample 
of  hawkmoths  for  study.  In  order  to  understand  the  general  characteristics  of  the  wings, 
Norris  performed  a  frequency  analysis  of  the  wing  using  a  scanning  laser  vibrometer 
(SLV).  Norris  vibrated  liberated  wings  from  the  hawkmoth  and  vibrated  them  at  various 
frequencies  using  an  SLV  to  measure  the  wing’s  response  in  both  air  and  in  vacuum.  The 
four  mode  shapes  he  determined  can  be  seen  below  in  Figure  1.6.  The  ratios  of  these 
frequencies  for  the  various  mode  shapes  show  the  relative  dynamic  stiffness  of  the  wing. 
These  mode  shapes  were  relatively  constant  across  a  large  sample  of  wings.  Norris 
identified  the  first  four  modes  as  flap,  feather,  saddle,  and  bisaddle  modes  respectively 
[6]. 


Figure  1.6:The  First  Four  Frequency  Modes  of  the  Hawkmoth's  Forewing. 


Norris  used  a  total  of  60  wings  to  perform  his  analysis.  While  the  mode  shapes 
remained  relatively  constant  from  specimen  to  specimen,  there  was  some  variation  [6]. 
Table  2  below  shows  the  results  of  Norris’  SLV  frequency  tests.  The  fact  that  there  is 
some  variation  in  the  results  is  significant  for  those  wishing  to  fabricate  a  bio-inspired 
wing  based  on  the  hawkmoth  because  it  means  that  the  wings  do  not  need  to  meet  exact 
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specifications  in  order  to  produce  an  artificial  wing  that  can  mimic  the  hawkmoth  wing. 
This  degree  of  tolerance  in  the  fabrication  of  such  an  MAV  is  important,  and  will  be 
shown  to  be  quite  significant  later  on. 

Table  1 .2:  Summary  of  the  Modal  Parameters  of  a  Hawkmoth's  Forewing  in  Air  and 
Vacuum  [6],  *  Based  on  10  samples  tested  in  near  vacuum 
**  Based  on  50  samples  tested  in  air 


Mode 

Structural 

Behavior 

Name 

Avg. 

Freq 

(Air)** 

[Hz] 

Avg. 

Freq 

(Vac)* 

[Hz] 

MR 

(Air)** 

[-1 

MR 

(Vac)* 

[-1 

Damp 

(Air)** 

[%1 

Damp 

(Vac)* 

[%1 

1 

SW 

Bending 

Flap 

60 

85 

1.0 

1.0 

5.0 

2.5 

2 

SW 

Torsion 

Feather 

84 

105 

1.4 

1.3 

5.0 

2.5 

3 

CW 

Bending 

Saddle 

107 

138 

1.8 

1.6 

5.0 

2.5 

4 

CW 

BiSaddle 

142 

170 

2.4 

2.2 

5.0 

2.5 

O’Hara  continued  on  the  work  of  Norris  in  hopes  of  developing  a  manufactured 
wing  that  could  be  used  in  an  MAV.  To  successfully  understand  the  structure,  it  is 
important  to  understand  the  geometry  and  material  proeprties.  O’Hara’s  work  looked  to 
provide  just  that.  For  his  research,  the  material  properties  of  the  wing  were  broken  down 
into  two  separate  parts,  the  vien  structure  and  the  membrane.  Both  the  membrane  and  the 
veins  were  then  characterized  using  a  variety  of  techniques.  Figure  1.7  below  shows  a 
venation  diagram  of  the  hawkmoth,  illustrating  the  names  and  locations  of  the  veins  in 
the  forewing  [11]. 
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Center  of  Pressure 


Costa  (C)  -  the  leading  edge  of  the  wing 

Subcosta(Sc)  -  second  longitudinal  vein  (behind  the  costa),  typically  unbranched 
Radius  (R)  -  third  longitudinal  vein,  one  five  branches  reach  the  wing  margin 
Media  (M)  -  fourth  longitudinal  vein,  one  to  four  branches  reach  the  wing  margin 
Cubitus  (Cu)  fifth  longitudinal  vein,  one  to  three  brances  reach  the  wing  margin 
Anal  Vein  (A)  -  unbranched  veins  behind  the  cubitus 


Figure  1.7:  A  Venation  Diagram  of  the  Manduca  Sexta  Forewing  [11]. 


In  order  to  determine  the  material  properties  of  the  veins,  it  was  important  to  first 
know  the  vein  dimensions.  Using  Computed  Tomography  (CT)  scan,  O’Hara  was  able  to 
develop  a  detailed  geometry  of  the  venation  pattern  of  the  hawkmoth,  as  well  as  take 
measurements  of  inner  and  outter  diameter  of  the  veins  [2].  Using  this  information,  he 
then  cut  out  individual  veins  from  the  liberated  wing.  Using  dynamic  forced  response, 
the  veins  were  then  tested  under  simple  cantilever  beam  conditions  to  predict  the  elastic 
modulus,  E,  of  the  structure.  This  was  done  using  laser  vibrometry  and  modal  analysis. 
A  section  of  the  vein  was  clamped  and  psuedo-randomly  vibrated.  The  frequency  of  the 
first  bend  was  measured  using  a  laser  vibrometer  (See  Figure  1.8  for  setup  and  results 
[2]). 
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Figure  1.8:  Experimental  Modal  Analysis:  Setup  and  Results  [2]. 


Using  a  simplified  finite  element  analysis  (FEA)  model  of  the  vein  and 
optimization  techniques,  a  value  for  E  of  the  FEA  model  was  optimized  so  that  the  the 
first  modal  response  frequency  observed  in  the  vibrometer  test  were  matched  in  the  FEA 
model.  This  methodology  has  been  used  previously  to  determine  various  composite 
material’s  properties  successfully  [12]. 

nmodes 

71=1 


60 f.n, 


- 1 


(l.i) 
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Equation  1.1,  depicts  the  simple  cost  function  used  to  tune  the  fundamental 
frequency  of  the  model  to  the  experiment.  By  iterating  on  the  elastic  modulus  variable  of 
the  FEA  model,  the  minimization  of  J  is  easily  realized  as  depicted  in  Equation  1.1. 
Using  this  method,  E  of  the  veins  was  detennined  to  be  6.8  GPa  [2]. 

O’Hara  used  a  Wyko  NT900  Optical  Profiler  that  was  able  to  determine  to  a  sub¬ 
nanometer  accuracy  of  the  shape  and  thickness  for  the  membrane  [2].  O’Hara  also 
implemented  instrumented  indentation  of  the  membrane,  a  common  practice  for 
detennining  the  mechanical  properties  of  thin  films  and  small  structural  features  [2]. 
Figure  1.9  describes  the  nano  indentation  process  used  by  O’Hara.  Using  twenty-five 
specimens  and  taking  multiple  measurements  for  each  specimen,  O’Hara  determined  the 
the  elastic  modulus  of  the  membrane  was  3.12  GPa  [2], 


Berkovich  Tip 


A  prescribed  load  is  applied  to  an  indenter  in 
contact  with  a  specimen.  As  the  load  is  applied, 
the  depth  of  penetration  is  measured.  The  area 
of  contact  at  full  load  is  determined  by  the 
depth  of  the  impression  and  the  known  angle  or 
radius  of  the  indenter.  The  hardness  is  found  by 
dividing  the  load  by  the  area  of  contact.  Shape 
of  the  unloading  curve  provides  a  measure  of 
elastic  modulus. 


P :  applied  load 

h :  indenter  displacement 

hr :  plastic  deformation  after  load  removal 

he  :  surface  displacement  at  the  contact  perimeter 


Figure  1.9:  The  Nanoindentation  Process  [2]. 


1.6.  Manufactured  Engineered  Wing  Materials 

The  development  of  FWMAV  requires  that  the  wings  be  able  to  produce 

significant  lift  with  the  least  amount  of  energy  input.  This  necessitates  a  wing  that  can 
serve  multiple  functions,  be  easily  controlled,  and  have  a  low  mass.  A  manufactured 
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wing  based  off  the  biological  wing  needs  to  be  able  to  match  the  physical  characteristics 
of  the  biological  wing,  as  well  as  be  repeatably  manufactured.  Keeping  this  in  mind, 
several  different  types  of  materials  were  looked  at  when  selecting  materials  to  be  used  in 
the  design  of  the  manufactured  MAV  wing.  Two  different  functions  for  the  materials 
were  examined,  use  in  the  vein  structure  and  use  in  the  membrane  of  the  manufactured 
wing.  Since  the  wing  would  be  primarily  experiencing  bending  in  flight,  the  goal  in 
selecting  a  material  was  that  the  flexural  stiffness  of  the  engineered  wing  would  match 
the  flexural  stiffness  of  the  bio  wing,  equation  1.2. 


E bio  1 bio  E on  n  I , 


eng1eng 


(1.2) 


Ultimately,  despite  looking  at  several  other  materials,  it  was  detennined  that  the  only 
material  that  is  readily  availible,  capable  of  achieving  this  task,  and  has  a  low  enough 
mass  to  be  used  as  the  vein  structure  for  the  engineered  wing  is  a  unidirectional  carbon 
fiber  composite.  Kapton  was  chosen  to  represent  the  membrane  based  on  its  strengh  and 
low  density. 


1.6.1.  Introduction  to  Unidirectional  Carbon  Fiber  Laminates 

When  performing  the  structural  analysis  of  the  engineered  wing,  it  is  of  the 
utmost  importance  to  understand  the  material  properties  of  the  unidirectional  carbon  fiber 
composite  that  make  up  the  vein  structure  of  the  model.  This  section  will  serve  to 
introduce  the  reader  to  the  basics  of  such  composites.  Unidirectional  carbon  fiber 
composites  are  the  most  common  example  of  continuous  fiber  composites,  composites  in 
which  the  fiber  is  continuous  throughout  the  length  of  the  entire  specimen. 
Unidirectional  carbon  fibers  are  continuous  fiber  composites  where  all  the  fibers  are 
aligned  in  one  direction,  hence  the  name  unidirectional. 
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The  nature  of  the  fibers  makes  such  unidirectional  carbon  fiber  composites 
anisotropic.  The  orientation  of  the  unidirectional  fibers  plays  a  significant  part  in  the 
material  properties  of  the  material.  The  orientation  of  a  fiber  is  the  angle,  9,  in  which  that 
fiber  forms  with  the  global  axis  coordinate  system.  Figure  1.10  serves  to  graphically 
define  the  lamina  coordinate  systems  (numbered)  and  the  global  laminate  coordinate 
system  (lettered). 

3,2 


Figure  1.10:  The  Principal  Directions  of  the  Lamina  (1,2,3),  and  the  Reference 
System  of  the  Laminate,  (x,y,z)  [13]. 


For  construction  purposes,  the  unidirectional  fibers  are  set  in  a  matrix,  usually 
epoxy,  and  fonned  into  thin  sheets  called  lamina,  or  plies  (See  Figure  1.11)  [14], 


Fibers 
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Laminates  are  made  by  stacking  the  lamina,  or  plies,  of  the  unidirectional  fiber  at 
various  fiber  orientations,  (See  Figure  1.12).  Laminates  are  beneficial  in  that  they  can  be 
arranged  in  a  near  infinite  number  of  orientations  and  ply  numbers  so  that  the 
composite’s  material  properties  can  be  tailored  to  specific  design  (See  Chapter  2). 


Figure  1.12:  Unidirectional  Lamina  Stacked  in  order  to  form  a  Laminate. 

1.6.2.  Unidirectional  High  Modulus  Thin  Ply  Laminates  as  MA  V  Wing  Structure 

Materials  such  as  steel,  titanium  and  plastic  polymers  were  found  to  be  either  too 
heavy  or  too  weak.  This  showed  the  necessity  for  a  material  with  a  high  specific 
modulus,  E/p.  Composites  composed  of  carbon  fiber  and  epoxy  resin  seemed  like  a 
natural  choice.  Not  only  do  they  possess  a  high  specific  modulus,  but  they  also  offer  the 
ability  to  tailor  the  number  of  plies  as  well  as  the  fiber  orientation  of  those  individual 
plies.  This  allows  for  more  control  over  how  the  material  will  respond  under  loading 
conditions,  allowing  for  a  fiber  orientation  and  ply  number  that  could  best  suit  the 
loading  conditions  applied.  Because  of  the  unidirectional  carbon-fiber  composite  is  not 
isotropic,  it  caused  some  significant  difficulties  when  performing  the  analysis.  The  issues 
associated  with  the  carbon  fiber  and  modeling  will  be  explained  in  subsequent  chapters. 
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Unidirectional  high  modulus  carbon  fiber  composites  have  already  proven  to  be 
useful  on  MAVS.  The  Harvard  Micro  Robotics  Research  Laboratory  has  used  a  pitch 
based  carbon  fiber  (XN-50A)  with  a  cyancate  ester  resin  (RS-3C)  produced  by  Tencate 
on  their  MAV  design  [15].  The  pitched  based  fibers  of  the  selected  composite  have  a 
high  strength  and  modulus  compared  to  other  types  of  available  composites  (See  Figure 
1.13  [16]).  This  makes  them  desirable  because  less  material  can  be  used  to  achieve  the 
same  strength  within  the  model,  but  at  a  lower  mass.  The  material  was  set  up  in  a 
[0/90/0]  orientation,  allowing  the  composite  to  have  stiffness  in  both  the  spanwise  and 
chordwise  directions  of  the  wing,  however  difficulties  were  found  in  the  manufacturing 
process  that  were  not  presented  in  the  literature.  These  issues  will  be  discussed  in 
subsequent  chapters. 


O  600 

s 


Tensile  Strength,  GPa 


Figure  1.13:  Composite  Fiber  Types 
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Since  the  carbon  fiber  used  by  Harvard  team  was  unavailable,  O’Hara  selected 
YSH-70A  fiber  for  use  as  the  vein  structure  in  the  engineered  wing,  material  properties 
shown  below  in  table  3. 


Table  1.3:  YSH-70A  High  Modulus  Fiber  Properties 


Fiber 

Tensile  Modulus 

Tensile  Strength 

Density 

Fiber  Dia 

Yield 

YSH-70A 

720  GPa 

3.6  GPa 

2.14  g/cm3 

7  pm 

125  g/ 1 000 m 

Previous  work  has  shown  that  the  inertial  properties  and  flexural  stiffness  of  the 
wing  play  an  important  part  in  the  dynamic  and  structural  response.  It  is  important  to 
understand  how  the  carbon  fiber  vein  structure  would  affect  an  engineered  wing.  Since 
the  carbon  fiber  material  itself  can  have  variations  in  the  ply  orientation  angle  and  the 
number  of  plies,  it  is  important  that  any  analysis  of  the  structure  consider  this.  It  is  also 
important  to  note  that  since  the  composite  is  anisotropic,  the  geometry  will  affect  the 
material  properties  at  different  points  along  the  structure. 

1.7.  Finite  Element  Approach 

Many  consider  FEA  to  be  one  of  the  most  important  structural  innovations  in 
recent  history.  It  allows  for  user  inputs  of  geometry  and  mass  to  perform  complex 
structural  analysis  much  quicker  than  would  be  possible  by  hand.  It  seems  only  natural 
that  the  finite  element  approach  should  be  used  in  the  analysis  of  an  engineered  MAV 
wing.  As  Travis  Sims  [17]  has  shown  with  his  work,  FEA  can  be  a  powerful  tool  to  solve 
a  problem,  but  it  can  also  lead  to  more  questions. 

1. 7.1.  A  Modal  Frequency  Approach 

Like  Norris,  Sims  understood  the  importance  of  the  structure  of  a  wing.  Sims 
objective  was  to  create  a  model  of  the  Manduca  Sexta’s  forewing  grounded  in 
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experimental  vibration  testing,  See  Figure  1.14  [17]  [7].  In  his  work,  as  with  all  FEA 
models,  there  are  two  competing  requirements:  maximizing  geometric  simplifications  to 
ensure  tractability,  and  minimize  unnecessary  deviations  from  the  physical  structure. 
Sim  developed  the  geometry  of  his  model  using  CT  imaging.  Since  the  real  material 
properties  of  the  bio  wing  were  unknown,  Sims’  generated  material  properties  so  that  his 
model  would  match  the  observed  modal  frequency  results  exhibited  by  a  liberated 
hawkmoth  wing  [17]. 


Figure  1.14:  Finite  Element  Model  of  Manduca  Sexta  Forewing  [17]. 

Sims’  model  yielded  similar  results  to  those  seen  by  Norris  in  his  modal  analysis 
except  for  the  third  mode,  See  Table  1.4  [17].  While  Sims  does  not  offer  a  clear 
explanation  for  this,  he  notes  that  the  mode  shapes  of  the  first  three  modes  examined  in 
the  model  were  the  same  as  for  the  biological  wing. 

Table  1.4:  Natural  Frequency  Results  Generated  by  Travis  Sims  [17]. 


Mode 

Experimental  Hz 

FE  Model  Hz 

Minimum  Difference 

1 

86  +/-  2 

84.6 

0.0% 

2 

106+/- 2 

106.1 

0.0% 

3 

155  +/-  2 

317.7 

102.4% 
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The  work  done  by  Sims  showed  that  material  properties  were  not  the  only 
important  aspect  in  matching  a  wing,  as  had  been  previously  suggested  by  Combes  and 
Daniels  [18].  Sims  came  up  with  the  following  areas  of  interest  that  should  be  looked  at 
when  determining  the  characteristics  of  a  bio  wing  [17]. 

1 .  The  material  properties. 

2.  The  geometry/shape  of  the  wings,  specifically  the  venation  pattern  and  vein 
structure. 

3.  The  natural  unique  camber  of  the  wing. 

In  an  effort  to  improve  his  model,  Sims  added  some  camber  to  his  wing  by 

applying  the  wing  outline  to  a  constant  camber  cylinder.  While  groundbreaking,  this 
camber  does  not  represent  the  actual  camber  found  on  the  wing.  However,  from  Sims’ 
results,  one  can  see  the  importance  that  camber  will  have  on  the  wing  (See  Figure  1.15). 


Figure  1.15:  Effect  of  camber  on  ton  for  the  Manduca  Sexta  Forewing  [17]. 

The  hawkmoth  is  not  the  only  insect  currently  being  analyzed  as  a  potential  for  a 
bio-inspired  MAV.  Marrocco,  Venkataraman  and  Demasi  have  investigated  the  use  of  a 
dragonfly’s  wing  and  have  developed  a  finite  element  model  [19].  For  their  model,  they 
assumed  a  planar  shape  and  used  material  properties  for  steel  to  represent  the  veins  and 
aluminum  for  the  membrane  since  material  properties  of  the  dragonfly  were  unavailable. 
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They  were  however  able  to  import  a  complex  geometry  of  the  wing,  and  vary  the 
thickness  of  the  wing  across  its  span  and  chord  length  (See  Figure  1.16). 


Figure  1.16:  FEA  Model  of  Marrocco  et  al  Dragonfly  Hindwing  [19]. 


Perfonning  a  modal  frequency  analysis,  Marrocco  et  al  were  able  to  see  the 
effects  the  veins  played  in  the  modal  shapes  of  the  wing.  A  second  set  of  runs  were  done 
to  compare  the  effects  of  mass  on  the  model  (See  Figure  1.17).  While  this  information 
provides  an  important  insight  into  the  dragonfly,  the  steel  and  aluminum  wing  is  not 
intended  to  mimic  the  dragonfly’s  wing,  or  an  engineered  wing  based  on  the  dragonfly. 
It  was  merely  to  investigate  the  mode  shapes  of  the  wing. 


Mode  ?  (Root  Con  drained  (. W)  f  Hv 


Figure  1.17:  Modal  Frequency  Analysis  of  Engineered  Dragonfly  Wing  [19]. 
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Further  literature  review  on  structural  analysis  of  FMAV  wings  shows  a 
significant  lack  of  results.  Analysis  on  manufactured  MAV  wings  is  still  in  its  infancy. 
One  recent  examples  that  could  be  found  included  work  done  by  Malik  and  Qureshi  [20]. 
While  the  work  proved  to  be  insightful,  the  structures  analyzed  did  not  have  the 
geometric  or  material  complexity  of  the  MAV  wing  that  is  the  subject  of  this  thesis. 
Malik  and  Qureshi  used  both  MATLAB  and  ANSYS  for  their  analysis.  Their  model  has  a 
solid  leading  edge  beam  with  a  trailing  membrane,  see  Figure  1.18.  This  figure  shows 
Malik  and  Qureshi ’s  model  in  ANSYS  (top)  and  MATLAB  (bottom)  undeformed  (left), 
and  the  first  modal  shape  [20].  The  frequency  for  the  first  mode  was  less  than  2  Hz.  It 
should  be  noted  that  this  is  significantly  lower  than  the  frequency  in  which  most  insects 
flap  [6]. 


Figure  1.18:  ANSYS  and  MATLAB  Models  of  Flapping  Wings  by  Malik  and  Qureshi 

[20]. 


22 


Sims’  areas  of  interest,  along  with  the  mass  properties  deemed  to  be  important  by 
O’Hara  and  Marrocco  et  al,  show  that  designing  a  MAV  wing  structure  is  not  an  easy 
task,  especially  considering  the  uncertainty  of  some  of  the  biological  wing  characteristics. 
It  is  easy  to  see  then,  how  a  better  understanding  of  an  engineered  wing  would  be 
beneficial.  An  analysis  of  O’Hara’s  engineered  MAV  wing  will  allow  for  better 
understanding  of  the  engineered  structure  and  how  it  varies  from  its  biological 
counterpart. 

1. 7.2.  A  Flexural  Stiffness  Approach 

Both  Sims  and  Norris  were  influenced  by  the  works  of  Combes  and  Daniels,  both 
biologists  at  the  University  of  Washington.  Combes  and  Daniel  observed  that  the  large- 
scale  deformations  observed  during  the  flight  of  insects  were  controlled  by  the 
architecture  of  the  wing  [18].  This  work  represents  one  of  the  few  works  that  showed 
investigations  into  the  structural  aspects  of  the  wings. 

The  parameter  that  Combes’  and  Daniels’  decided  was  important  to  investigate 
was  flexural  stiffness.  They  define  it  as  “the  composite  measure  of  the  overall  bending 
stiffness  of  a  wing;  it  is  the  product  of  the  material  stiffness  (E,  which  describes  the 
stiffness  of  the  wing  material  itself)  and  the  second  moment  of  area  (I,  which  described 
the  stiffness  generated  by  the  cross  sectional  geometry  of  the  wing)”  [18].  The 
mathematical  equation  for  a  beam  is  shown  below  in  Equation  1.3. 

FT  FL' 3 

EI  =  IT  <U) 

Here  the  parameter  [L]  is  the  effective  beam  length,  [5]  is  the  wing  displacement 
at  the  given  position  of  force  application,  and  [F]  is  the  applied  force.  For  their  case, 
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Combes  and  Daniel  used  the  following  equation,  Equation  1.4,  to  solve  for  the  second 
moment  of  area.  In  this  case  [vv]  is  the  width,  and  [/]  is  the  thickness  of  the  wing. 


(1.4) 
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Combes  and  Daniel  set  up  a  way  to  test  this  parameter  based  on  experimental  and 
FEA  results  (like  Sims)  which  is  represented  in  Figure  1.19.  The  figure  shows  the  wing 
displacement  tests  used  to  measure  the  flexural  stiffness  of  the  wing.  Figure  1 .20  shows 
the  finite  element  model  representative  of  Combes’  and  Daniels’  flexural  stiffness 
experiments.  A  similar  method  will  be  employed  in  this  thesis  for  analysis  of  the 
engineered  wing. 


A 


70%  span 


70%  chord 


Figure  1.19:  Combes’  and  Daniel's  Initial  Flexural  Stiffness  Investigations  [18]. 
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A 


Vein  material  stiffness  (N  nr2) 


Figure  1.20:  The  FEA  Model  with  Results  Developed  by  Combes  and  Daniels 


Through  their  investigation,  Combes’  and  Daniels’  were  able  to  measure  the 
stiffness  of  the  wing,  and  it  is  from  this  experiment  that  they  were  able  to  measure  the 
chordwise  and  spanwise  stiffness  of  several  different  wings.  Looking  to  standardize  the 
size  of  the  specimens  tested,  several  different  species  wings  were  tested.  Plotting  the 
spanwise  and  chordwise  flexural  stiffness  versus  the  wing  span  and  chord  length 
respectively  for  these  tests,  Combes’  and  Daniels  came  up  with  the  conclusion  that  size 
scaling  was  the  dominant  factor  in  determining  overall  flexural  stiffness  [18].  It  is 
important  then,  that  when  manufacturing  a  FMAV  for  the  wing  to  follow  similar 
spanwise  and  chordwise  characteristics  to  those  found  by  Combes’  and  Daniels’,  Figure 
1.21. 
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Figure  1.21:  Results  from  Combes’  and  Daniels’  Measurement  of  Flexural  Stiffness  [18]. 


Combes’  and  Daniels’  approach  to  analyzing  wings  has  shown  some  common 
characteristics  are  shared  amongst  different  species.  Their  FEA  model  like  Sims  original, 
did  not  characterize  the  camber  of  the  wing.  While  flexural  stiffness  measurements 
should  provide  yet  another  tool  to  compare  the  FEA  model  against  the  engineered  wing,  a 
more  complex  model  than  is  needed  to  handle  the  frequency  modal  analysis  as  well  as 
point  loads  measuring  flexural  stiffness.  This  is  important  to  consider  when  developing 
an  FEA  model  for  the  engineered  wing. 


1.8.  Objective  and  Document  Ovierview 

Looking  at  what  work  has  been  done,  it  can  be  seen  that  the  area  of  MAV 

research  still  has  a  long  way  to  go.  The  objective  of  this  research  is  to  perform  a 
structural  analysis  on  an  engineered  Manduca  Sexta  forewing  with  a  composite  vein 
structure.  The  composite  material  was  examined  to  determine  issues  associated  with  the 
use  of  an  anisotropic  material  in  MAV  wing  contraction.  Also,  as  part  of  the  analysis  an 
FEA  model  was  created  that  had  the  ability  to  do  the  following: 
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1 .  Accurately  predict  mode  shapes  and  frequencies  of  a  modal  frequency  analysis. 

2.  Accurately  model  displacements  due  to  point  loads. 

3.  Be  easily  modified  to  incorporate  any  future  design  changes. 

The  following  chapters  will  discuss  how  this  process  was  carried  out,  as  well  as 
any  issues  encountered  and  possible  solutions.  Chapter  II  will  detail  some  of  the  theory 
used  in  this  project  in  dealing  with  the  carbon  fiber  composite  material,  and  the  use  of 
finite  element  method.  Chapter  III  will  explain  the  construction  of  the  composite 
material  and  the  engineered  wing,  and  go  over  the  analysis  of  the  composite  material. 
Chapter  IV  will  go  over  how  the  FEA  model  of  the  engineered  wing  was  developed,  as 
well  as  experimentation  done  to  validate  the  model.  Chapter  V  will  display  the  results  of 
the  FEA,  as  well  as  provide  a  discussion  on  those  results.  Finally,  Chapter  VI  will 
provide  conclusions  and  a  summary  of  the  work  done  in  this  thesis. 
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II.  Carbon  Fiber  Composite  Construction  and  Theory 


This  chapter  serves  to  demonstrate  an  understanding  of  the  mathematics  involved 
in  the  analysis  performed  for  this  project.  The  two  main  areas  of  focus  include  the 
material  properties  of  carbon  fiber  composite  material,  which  forms  the  vein  structure  of 
the  wing  (See  Figure  2.1),  and  an  understanding  of  FEA  frequency  analysis,  and  how  the 
composite  was  modeled  so  that  such  an  analysis  could  be  performed.  Figure  2.1  is  a 
diagram  depicting  the  wing  structure  analyzed  and  manufactured.  The  reader  should 
notice  how  the  individual  vein  shapes  cut  through  the  associated  composite  ply 
orientations.  It  is  important  to  understand  how  this  was  handled,  and  this  chapter  will  go 
over  some  of  the  techniques  used. 


Figure  2.1:  The  Manufactured  Wing  with  Lines  Representing  the  0°  and  90°  Layout  of 

the  Carbon  Fibers. 

2.1.  Unidirectional  Carbon  Fiber  Material  Construction 

This  section  will  serve  to  define  the  material  properties  of  the  YSH-70-A/RS-3C 

carbon  fiber  composite  used  in  the  vein  structure  of  the  wing,  show  how  the  material  was 
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constructed,  as  well  as  go  over  the  theory  behind  calculations  that  were  used  in  the 
analysis  of  the  composite  itself  and  in  the  analysis  of  the  manufactured  wing. 


The  lamina  material  properties  are  provided  by  the  manufacturer  as  well  as 
confirmed  by  experimental  methods,  however  the  laminate  needed  to  be  constructed  and 
its  properties  calculated.  The  laminate  properties  varied  from  the  lamina  properties  due 
the  variation  of  the  ply  orientation  of  the  libers.  In  the  case  of  the  engineered  wing,  a  3- 
ply  lamina  was  arranged  in  a  [0/90/0]  orientation,  as  shown  in  Figure  2.2. 


Figure  2.2:Diagram  of  the  [0/90/0]  Composite  Laminate. 


Two  main  sets  of  calculations  were  performed  during  the  analysis  of  the 
composite,  the  lamina  transformation  equations,  which  account  for  the  change  in  material 
properties  of  the  lamina  based  upon  the  fiber  orientation  angle,  and  the  Halpin-Tsai 
equations,  which  account  for  changes  in  the  material  properties  due  to  short  fibers.  Short 
fibers  occur  when  the  length  of  the  fiber  is  approximately  less  than  100  times  the 
diameter  of  the  fiber,  and  have  a  force  component  acting  along  the  length  of  the  fiber. 
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Such  fibers  occur  in  the  manufactured  wing  due  to  the  curvature  of  the  vein  structure  cut 
from  the  single  sheet  of  the  [0/90/0]  composite  as  shown  in  Figure  2.3. 

2.1.1.  Laminate  and  Wing  Construction  at  AF IT 

Part  of  the  analysis  was  detennining  how  it  was  constructed.  The  laminate  used 
in  the  fabrication  of  the  engineered  wing  was  constructed  at  AFIT.  This  was  specifically 
done  from  an  economical  efficiency  point  of  view.  The  pre-preg,  lamina  where  the  fibers 
are  already  impregnated  in  a  partially  cured  matrix,  was  used  by  several  groups  of 
students,  all  constructing  small  samples  of  the  composite  at  various  times,  in  various  ply 
orientations  and  various  numbers  of  plies.  It  was  deemed  impractical  to  have  every 
student  group  purchase  small  samples  of  cured  carbon  fiber  for  every  application.  Since 
3 -ply  laminate  is  not  common,  this  kept  different  groups  of  students  from  having  to 
special  order  their  composite  which  can  be  very  expensive  for  the  small  amount  of 
carbon-fiber  used.  It  also  had  the  advantage  of  allowing  AFIT  students  to  lay-up  and 
cure  their  laminate  within  a  day  as  opposed  to  having  to  order,  and  wait  an  unknown 
amount  of  time  for  it  to  be  shipped. 

The  composite  arrives  at  AFIT  in  rolls  of  pre-preg.  Since  the  pre-preg  must  be 
stored  below  32°F  in  order  to  prevent  the  matrix  from  curing  prematurely,  the  roll  of  pre- 
preg  is  cut  with  a  straight  edge  and  razor  blade  into  8”  x  1 1”  sheets  (See  Figure  2.3)  so 
that  they  are  able  to  be  stored  in  a  common  household  freezer  at  AFIT  until  needed. 
These  sheets  will  make  up  the  individual  lamina  of  the  composite.  Note:  The  roll 
pictured  on  top  in  Figure  2.3  is  not  the  actual  carbon  fiber  pre-preg  used  at  AFIT  because 
the  composite  has  already  been  cut  up.  The  actual  pre-preg  is  pictured  on  the  bottom,  and 
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is  sandwiched  between  thin  sheets  of  wax  paper  on  one  side  and  a  thin  film  of  plastic  on 
the  other  side. 


Figure  2.3:  A  Diagram  of  the  Roll  of  Pre-preg  Carbon  Fiber  Composite  (Top),  and  Actual 
Composite  Sheet  with  Protective  Film  (Bottom). 

Once  the  plies  are  ready  to  be  stacked,  the  protective  plastic  film  and  wax  paper 
sheets  are  removed  from  each  ply.  These  plies,  with  the  help  of  a  straight  edge,  are 
stacked  one  on  top  of  another  in  the  desired  [0/90/0]  orientation. 

This  new  uncured  laminate  is  then  prepared  to  be  placed  in  the  heat  press.  The 
uncured  laminate  is  placed  in  between  Teflon,  metal  plates,  and  cardboard  plates  as 
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shown  in  Figure  2.4.  The  heat  transfer  paper  and  particleboard  used  are  to  ensure  that  the 
force  from  the  heat  press  is  evenly  distributed  to  the  composite,  and  the  steel  and  metal 
plates  allow  for  the  composite  to  remain  flat.  The  non-porous  TEFLON  ensures  that  the 
composite  does  not  adhere  to  the  metal  plates  during  or  after  the  curing  process,  and  the 
porous  TEFLON  allows  for  excess  matrix  material  to  bleed  out  without  adhering  to 
anything. 


1. 

Heat  Transfer  Paper, 

2. 

Steel  Heat  Press  Plates 

3. 

Particle  Board 

4. 

Metal  Plates 

5. 

Non-Porous  TEFLON 

6. 

Porous  Teflon 

7. 

Carbon  Fiber  Composite 

Figure  2.4:  Set-up  of  the  Material  for  the  Heat  Press. 

The  pre-preg  stack  is  then  placed  into  a  LPKF  Multipress  S,  Figure  2.5,  a  heat 
press  that  is  capable  of  providing  both  pressure  and  heat  required  to  cure  the  pre-preg 
material.  The  material  was  pressed  initially  for  10  minutes  at  30°C  and  30N/cm  ,  then 
for  120  minutes  at  192°C  and  100  N/cm".  Pressure  was  removed  and  the  composite  was 
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allowed  to  cool  at  room  temperature  for  several  hours  before  being  removed  from  the 
oven,  Figure  2.5.  This  process  cures  the  laminate  so  that  it  is  then  ready  to  be  used. 


Figure  2.5:  LPKF  MultiPress  S. 


Once  the  carbon  fiber  had  been  cured,  the  composite  is  ready  to  be  cut  into  either 
the  wing  or  other  test  specimens.  Due  to  the  small  nature  of  the  veins  or  the  specimens,  a 
very  precise  method  is  required  for  manufacturing.  This  was  done  using  a  laser  cutter, 
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and  was  perfonned  at  several  locations,  the  Sensors  Directorate  at  WPAFB,  or  Mount 
Laser  and  Photonics  Center  (MLPC),  in  Miamisburg,  OH,  Figure  2.7.  While  the 
locations  and  the  lasers  used  for  the  cutting  were  different,  the  general  process  is  the 
same. 

The  single  sheet  of  carbon  fiber  is  placed  into  the  laser  cutting  area.  A  straight 
edge  is  used  as  a  guide  to  line  up  the  carbon  fiber  sheet  within  the  laser  itself  by  ensuring 
that  one  edge  of  the  sheet  is  placed  up  against  the  straight  edge.  The  carbon  is  held  in 
place  by  a  vacuum  being  pulled  from  below  the  flat  cutting  surface  of  the  laser.  A  .dxf 
file  is  created  previously  in  a  computer  automated  design  (CAD)  program  that  represents 
the  path  that  the  laser  will  follow  to  cut.  This  file  is  imported  into  the  computer  that 
controls  the  laser.  The  settings  of  the  laser  should  be  arranged  so  that  the  laser  is  able  to 
cut  through  the  carbon  fiber,  but  also  so  that  the  carbon  fiber  does  not  become  burnt 
during  the  process.  This  is  important  because  as  the  composite  burns,  the  integrity  of  the 
material  is  compromised.  Once  these  settings  have  been  detennined,  the  composite  is 
then  cut.  For  the  case  of  this  project,  several  sample  parts  were  made  and  examined  by 
MLPC  in  order  to  ensure  the  effectiveness  of  the  cut,  and  the  effects  of  the  cut  on  the 
composite  material.  This  ensured  that  as  the  material  was  cut  with  the  laser,  no  adverse 
affects  to  the  material  occurred,  such  as  overly  heated  edges,  which  could  cause  the 
material  to  bum  and  degrade  the  material  properties. 
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a)  LPKF  Laser  at  WPAFB  b)  Laser  at  MLPC 

Figure  2.7:  Laser  Cutters  used  a)  LPKF  PhotoLaser  S  at  WFAFB,  and  b)  In-House  Laser 

at  MLPC. 

Finally  the  Kapton  membrane  is  applied  to  the  wing.  This  is  done  using  3M  45 
spray  on  adhesive  that  is  sprayed  into  a  small  container,  which  is  then  applied  via  a  paint 
brush  by  hand  to  each  of  the  carbon  fiber  veins  in  the  wing.  The  membrane  is  then 
placed  onto  the  vein  structure,  and  the  adhesive  is  allowed  to  dry.  Once  the  adhesive 
dries,  excess  membrane  is  trimmed  away,  and  the  wing  is  then  ready  for  use,  Figure  2.8. 
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Application  of  the  Adhesive 


Membrane  Application 


Excess  Membrane  Removal 


Figure  2.8:  Application  of  the  Kapton  Membrane. 


Vein  Structure  Membrane 


Figure  2.9:  Completed  FMAV  Wing. 


2.2.  YSH-70-A/RS-3C 

This  section  will  serve  to  establish  the  baseline  material  properties  of  the  YSH- 
70-A/RS-3C  carbon  fiber  composite  used  in  the  vein  structure  of  the  wing.  The  baseline 
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material  properties  were  provided  by  the  manufacturer,  Appendix  A,  and  had  been 
confirmed  by  O’Hara,  who  had  previously  tested  20-ply  specimens  of  the  composite 
using  ASTM  D  3039  and  ASTM  D  3518  standards.  These  tests  were  done  in 
coordination  with  the  University  of  Dayton  Research  Institute  and  Air  Force  Research 
Laboratory  Materials  Directorate  on  ply  orientations  of  0,  90  and  +/-  45  degrees  test 
specimens  [2],  The  tests  would  allow  for  the  calculations  of  the  moduli  in  the  principal 
and  secondary  axes,  as  well  as  the  shear  direction.  Table  2.1  shows  the  resulting  lamina 
material  properties  as  determined  by  the  tests,  and  is  compared  with  the  properties  given 
by  the  manufacturer  [2], 

Table  2. 1 :  Material  Properties  developed  by  O'Hara  and  Manufacturer  [2]. 


YSH-70-A/RS-3C  Lamina  Material  Properties 

20  Plv 

Manufacturer 

%  Difference 

Ei 

4.15  E+ll  Pa 

4.20  E+ll  Pa 

1.20 

e2 

5.52  E+09  Pa 

5.51  E+09  Pa 

0.18 

Vl2 

3.00  E-01 

2.80  E-01 

7.14 

Gl2 

4.85  E+09  Pa 

4.83  E+09  Pa 

0.41 

While  there  is  a  slight  difference  in  material  properties  based  on  the  20  ply 
compared  to  the  manufacturers  specifications,  such  a  small  error  is  common  between  test 
specimens  and  the  material  properties.  Figure  2.10  [1]  shows  the  stress-strain  curve  for 
the  5  samples  of  the  20-ply  test  perfonned  by  O’Hara.  These  results  were  averaged  to 
arrive  at  the  given  modulus. 
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Figure  2.10:  20-Ply  Tensile  Test  Stress-Strain  Curve 


The  results  from  these  tests  confirm  that  the  base  material  properties  provided  by 
the  manufacturer  were  mostly  correct.  Therefore,  it  was  the  manufacterer’s  properties 
listed  in  Table  2. 1  above  that  were  used  in  further  analysis. 


2.3.  Lamina  Engineering  Constants 

In  order  to  perform  the  transfonnation  analysis  required  to  follow  the  veins 
curvature  as  previously  mentioned,  an  angular  reorientation  had  to  be  created  using  the 
geometry  shown  in  Figure  2.9.  These  transformations  account  for  the  fact  that  the 
laminate  was  constructed  in  a  [0/90/0]  orientation,  and  are  used  to  solve  for  the  material 
properties  of  the  lamina  for  instances  when  the  segments  considered  are  not  orientated  to 
0°. 
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Figure  2.11:  Off  Axis  Lamina  Under  Tensile  Stress,  with  Global  (x,y)  and  Material  (1,2) 

Axes 


The  material  properties  in  the  local  coordinate  system  (x,y,z)  are  as  follows  for  the  case 
of  a  beam: 


Ei 

m4+n4|^ 

e2\ 


(2.1) 


Where  Ex  is  the  elastic  modulus  in  the  x  axis,  Ei  is  the  elastic  modulus  of  the 
composite  at  0°,  E2  is  the  modulus  of  the  composite  at  90°  orientation,  m  is  the  cosine  of 
the  ply  orientation,  and  n  is  the  sine  of  the  ply  orientation.  It  is  important  to  recognize  the 
effect  that  the  ply  angle  has  on  this  equation. 

To  illustrate  the  effects  that  the  ply  orientation  has  on  material  properties  of  the 
lamina,  the  Modulus  in  Tension  is  plotted  against  the  lamina  orientation  for  YSH- 
70A/RS-3C,  (See  Figure  2.12). 
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Figure  2.12:  The  Effects  of  Tension  Modulus  vs.  Lamina  Orientation  Angle  of  the  Fiber 

for  YSH-70-A/RS-3C. 


It  can  be  seen  that  there  is  a  drastic  decline  in  modulus  of  the  material  as  the  ply 
orientation  deviates  from  0°.  This  is  of  extreme  importance,  especially  for  the  use  of  a 
thin  3-ply  laminate.  Considerable  amounts  of  time  were  spent  in  understanding  the 
effects  of  this  while  performing  the  analysis  as  will  be  shown  later  in  this  chapter. 


2.4.  Halpin-Tsai  (Short  Fibers) 

Since  the  vein  structure  being  fabricated  for  the  engineered  wing  is  very  narrow 
and  cut  across  the  individual  fibers,  short  fibers  that  may  have  occurred  were  accounted 
for  by  use  of  the  Halpin-Tsai  Equations,  Equation  2.5  [21].  Short  fibers  are  defined  as 
fibers  where  the  length  of  the  fiber  is  not  significantly  larger  than  the  diameter  of  the 
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fiber.  While  there  is  no  defined  length  to  diameter  ratio  for  a  fiber  to  be  defined  as 


short”,  the  accepted  ratio  is  typically  around  100  [21]. 


P  1  +  Qrjcp 


(2.5) 


Pm  1  -  C  V<t> 


The  engineering  constants,  Ef  and  Em,  were  directly  substituted  for  Pm  and  Pf  in 
Equation  2.5.  Where  Ef  is  the  modulus  of  the  fiber,  Em  is  the  modulus  of  the  matrix, 
obtained  from  the  manufacturer,  Appendix  A.  For  the  case  of  circular  fibers,  which  were 
present  in  the  wing,  the  following  factors  were  inserted  into  Equation  2.5.  Equations  2.6 
and  2.7  show  these  factors,  where  1  is  the  length  of  the  fiber  in  the  1 -direction,  t  is  the 
thickness  of  the  tape  and  w  is  the  width  of  the  fiber  in  the  2-direction. 


Ceu  =  2 m 

<e22  =  2{w/t) 

Con  =  1 


(2.6) 


(2.7) 


For  the  case  of  circular  fibers,  w=t.  It  should  be  noted  that  for  cases  where  the 
length  of  fiber  is  sufficiently  long,  no  significant  difference  in  material  properties  is 
occurred  through  the  use  of  the  Halpin-Tsai  equation.  Using  these  equations,  the  material 
properties  were  solved  for  the  lamina. 
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2.5.  Finite  Element  Analysis 

A  significant  amount  of  the  analysis  of  the  engineered  wing  was  perfonned  using 
FEA.  An  FEA  model  was  constructed  by  developing  a  MATLAB  program  that  is 
capable  of  generating  an  input  file  to  be  submitted  to  ABAQUS,  a  finite  element  model 
solver.  This  technique  was  used  in  order  to  more  easily  capture  the  complex  nature  of  the 
wing.  Frequency  modal  analysis  of  the  wing  was  performed,  as  well  as  point  load  tests  on 
the  wing.  The  following  sections  go  through  the  theory  and  processes  required  to 
perfonn  the  analysis  and  experimentation. 

2.5.1.  Finite  Element  Frequency  Analysis  Theory 

When  performing  the  frequency  analysis  of  the  wing,  it  is  important  to  understand 
the  fundamental  concept  behind  what  is  being  done  within  the  finite  element  analysis 
program.  While  the  basic  theory  behind  finite  element  analysis  when  dealing  with  static 
loading  is  readily  available,  the  theory  behind  using  FEA  for  frequency  analysis  is  not  as 
readily  available.  Therefore,  only  the  basics  behind  the  frequency  analysis  will  be 
covered  in  this  section. 

The  modal  frequencies  of  structure  are  based  on  the  structures  stiffness  and  mass 
properties.  In  the  finite  element  model,  the  mass  and  stiffness  are  represented  at  a  finite 
set  of  second  order  differential  equations  in  the  time  domain,  [t]  (equations  of  motion, 
Equation  2.8)  as  follows: 

K  ■  U  +  M  •  d2U/dt2  =  0  <2'8> 

Where  [K]  represents  the  stiffness  matrix  of  the  finite  element  model,  [U]  represents  the 
vector  of  nodal  displacements,  with  six  degrees  of  freedom  (three  rotation,  and  three 
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translation).  [M]  represents  the  finite  element  mass  matrix.  The  nodal  displacement 
vector  is  set  to  the  equation  below: 

(2.9) 

U  =  OeIC0t 

Substituting  the  equation  of  motion  (Equation  2.10),  obtains  the  following  eigenvalue 
problem. 

[K  -  co2-M]<D  =  0  <2-10> 

['!’]  is  the  resulting  mode  shape  (eigenvector),  and  [to]  is  the  associated  frequency 
(eigenvalue).  ABAQUS  was  used  to  calculate  and  display  the  mode  shape  and  associated 
frequency.  It  is  important  to  note  that  the  associated  modal  shape  and  frequency  is  also 
associated  with  the  boundary  condition,  which  may  not  be  evidently  clear  from  the 
equations  above.  The  boundary  condition  for  this  project  was  clamped,  prohibiting 
displacement  in  all  six  degrees  of  freedom.  Mathematically,  this  would  be  evident  in  the 
vector  of  nodal  displacement,  [U],  where  the  values  for  the  clamped  nodes  would  be  set 
at  zero. 

2.5.2.  Effective  Moment  of  Inertia 

Since  beam  elements  were  used  in  this  project,  properly  modeling  the  nature  of 
the  composite  is  important  since  the  material  does  not  have  constant  material  properties 
through  the  entire  cross  section.  The  effective  moment  of  inertia  allows  for  one  material 
property  to  represent  the  element  while  changing  the  dimensions  of  the  element’s  beam 
profile  to  maintain  its  physical  characteristics.  In  order  to  capture  this,  the  effective 
modulus  of  the  beam  cross  section  is  taken  for  each  beam  element  used  to  model  the 
composite  vein  in  the  model. 
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For  the  given  section  of  the  vein,  the  average  angle  is  calculated  and  the  ply 
orientations  for  each  element  are  detennined.  Using  the  code  in  Appendix  B,  the  material 
properties  for  each  lamina  are  detennined.  Since  only  one  set  of  material  properties  from 
the  beam  needs  to  be  used,  the  material  properties  from  the  lamina  with  the  highest  axial 
modulus  was  chosen  to  represent  the  entire  beam  with  the  exception  of  the  density,  which 
will  be  explained  in  Section  2.5.3.  Using  the  material  properties  from  the  ply  with  the 
highest  axial  modulus  for  the  element  allows  for  maximum  width  of  the  beam  profile  in 
the  model  to  be  viewed  within  the  model  as  the  same  physical  width  of  the  carbon-fiber 
at  the  same  point  on  physical  specimen.  Since  there  is  a  large  difference  in  the  stiffness 
of  the  material  based  on  fiber  direction,  this  prevents  the  maximum  width  of  the  beam 
element  profile  from  becoming  too  large  so  that  when  the  beam  profile  is  displayed  in 
ABAQUS  one  can  easily  see  a  shape  similar  to  the  actual  wing  and  not  a  vastly  distorted 
image. 

A  ratio  of  the  axial  moduli  of  the  other  two  plies  to  the  maximum  modulus 
determines  the  reduction  of  area  needed  by  those  plies  to  retain  the  physical 
characteristics  of  the  beam.  Equation  2.11  shows  the  equation  used  to  perfonn  the 
transfonnation. 

Nx= ^  (2.1.) 

xjnax 

Where  Expiy  is  the  elastic  modulus  of  the  lamina  in  the  local  x  direction.  Ex_max  is  the 
maximum  elastic  modulus  of  the  lamina  in  the  local  x  direction  for  the  laminate.  Nx  is 
the  axial  modulus  ratio  of  the  xth  ply. 
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The  axial  modulus  ratio  is  used  to  change  the  area  of  the  ply  so  that  when  the 
singular  material  property  is  applied,  the  beam  section  maintains  the  same  axial  modulus. 
In  order  to  keep  the  effective  moment  of  inertia  the  same  for  the  element,  the  distance  of 
the  ply  from  the  neutral  axis  needs  to  be  kept  constant,  therefore  when  changing  the  area 
of  the  height  and  location  of  the  ply  was  held  constant  and  only  the  width  was  modified. 
The  effective  width  of  each  of  the  two  other  plies  was  simply  the  width  of  the  beam 
section  multiplied  by  their  individual  axial  modulus  ratios,  Equation  2.12. 

^Kiew  Nx  %  Worig  (2.12) 

Where  Wnew  piy  is  the  new  effective  width  for  the  lamina  in  the  model’s  cross  section,  and 
Worig  is  the  original  cross  section  width  of  the  lamina  in  the  model. 

Implementing  the  effective  moment  of  inertia  calculation  changes  the  beam 
elements  profile  from  a  rectangular  shape  to  a  shape  similar  to  an  I-beam  for  most 
elements.  Figure  2.13  shows  a  representation  of  a  section  of  the  three  ply  carbon-fiber 
material  and  the  new  beam  section  that  would  result  based  on  the  effective  moment  of 
inertia  calculations. 
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Figure  2.13:  Diagram  Showing  the  Effect  on  the  Beam  Cross  Section  after  the  Effective 
Moment  of  Inertia  Calculations  is  Applied. 

2.5.3.  Effective  Density 

While  the  stiffness  of  the  element  represents  the  physical  characteristics  of  the 
carbon  fiber,  the  mass  of  the  element  is  changed  due  to  the  change  in  area  if  the  density  is 
kept  the  same.  In  order  to  account  for  this,  the  density  of  the  material  for  the  element  is 
changed  so  that  the  mass  of  the  element  is  the  same  as  it  was  prior  to  the  change  in  area 
taking  place.  This  is  a  simple  calculation  based  on  the  ratios  of  the  cross  sectional  area  of 
the  new  beam  section  compared  to  the  physical  cross  sectional  area  of  the  wing  at  that 
point,  Equation  2.13 


^  ^  ^  A0rig 

Pnew  Porig  x  ~A 

A new 


(2.13) 
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Where  pnew  is  the  new  density  for  the  element,  porig  is  the  original  density  of  the 
composite,  Aorig  is  the  average  physical  cross  section  of  the  manufactured  wing  at  the 
same  location  as  the  element,  and  Anew  is  the  new  cross  section  of  the  beam  element  after 
the  effect  moment  of  Inertia  calculations  have  been  performed. 
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III.  Carbon  Fiber  Composite  Analysis  and  Results 


This  chapter  will  serve  to  describe  the  work  done  to  analyze  the  material 
properties  of  the  carbon  fiber  composite  used  for  the  vein  structure  of  the  engineered 
wing.  Due  to  the  requirement  for  a  low  mass,  the  composite  was  chosen  to  serve  as  the 
vein  structure  due  to  its  high  strength  and  low  density.  Previous  work  done  at  Harvard 
showed  that  the  composite  was  a  viable  alternative  to  metals  or  polymers,  and  suggested 
that  a  [0/90/0]  orientation  was  useful  in  that  it  gave  the  composite  the  ability  to  be  stiff  in 
the  spanwise  and  chordwise  direction  [15].  However,  some  of  the  issues  encountered  in 
this  project  were  not  documented  in  the  literature.  It  was  important  to  understand  how 
these  issues  affected  the  material  properties  of  the  composite,  and  ultimately  in  the 
engineered  wings  themselves. 

One  can  observe  that,  as  shown  in  Figure  3.1,  the  actual  orientation  of  the  fiber  in 
a  particular  vein  is  predicated  upon  the  veins  curvature  at  a  coordinated  point.  This  make 
the  through  thickness  material  properties  per  segment  anisotropic.  Thus,  it  became 
necessary  to  not  only  develop  a  method  to  handle  this  orientation  requirement,  but  to 
make  sure  that  any  material  property  was  correctly  defined.  This  chapter  serves  to 
evaluate  the  unidirectional  carbon  fiber  composite  used  in  the  manufacturing  of  the 
engineered  MAV  wing.  The  vein  structure  of  the  wing  was  cut  out  of  a  single  sheet  of 
composite  material  made  of  three  plies  of  unidirectional  carbon  fiber  stacked  in  a  [0/90/0] 
orientation  as  shown,  Figure  3.1.  Evidence  suggested  some  issues  with  the  carbon  fiber 
composite  material  properties,  and  an  effort  was  made  to  understand  and  eventually 
characterize  material  variations  before  beginning  the  analysis  of  the  wing. 
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Local  Axis 

Figure  3.1:  The  Manufactured  Wing  with  Lines  Representing  the  0°  and  90°  Layout  of 

the  Carbon  Fibers. 


3.1.  Laminate  Material  Modal  Frequency  Experimentation 

Prior  to  the  analysis  of  the  engineered  wing,  there  were  several  questions  raised 

about  the  validity  of  the  unidirectional  carbon  fiber  composite  material  properties  used  in 
the  creation  of  the  engineered  MAV  wing  stemming  from  the  results  of  the  previous 
tensile  tests,  and  the  lack  of  literature  on  the  use  of  3-ply  specimens.  Since  the  carbon 
fiber  composite  is  used  to  construct  the  veins  of  the  engineered  wing,  Figure  3.1,  and 
provide  most  of  the  structure’s  mass  and  stiffness,  it  was  important  to  examine  any  of 
these  concerns.  Therefore,  an  analysis  of  the  carbon  fiber  composite  was  performed  in 
conjunction  with  the  analysis  performed  on  the  engineered  wing. 

3.1.1.  Laminate  Material  Sample  Testing 

Wishing  to  use  a  non-destructive  method,  a  modal  frequency  test  was  perfonned 
on  carbon  fiber  composite  specimens  that  were  40  mm  long  by  5  mm  wide  in  order  to 
examine  the  material  properties  and  determine  the  validity  of  the  given  material 
properties,  since  only  three  plies  would  be  used  to  construct  the  laminate.  These 
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specimens  were  three-ply  laminates  with  fibers  at  a  [0/90/0]  orientation.  The  specimens 
are  significantly  smaller  than  the  120  mm  by  25  mm  specimens  tested  by  O’Hara  when 
selecting  the  carbon  fiber  composite  for  the  wing  [2],  and  closer  to  the  size  of  the  veins 
that  are  present  on  the  engineered  wing. 

For  the  modal  frequency  test,  the  specimens  were  clamped  at  the  base,  and 
vibrated  with  airborne  excitation  in  a  pseudo  random  manner,  Figure  3.2.  This  modal 
frequency  test  is  the  same  method  that  was  used  to  determine  the  characteristcs  of  the 
vein  material  properties,  and  was  also  used  by  DeLeon  (See  Section  1.5,  and  Equation  2 
for  more  details)  [1]. 


Figure  3.2:  40mm  x  5  mm  Beam  with  Airborne  Excitation. 
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In  order  to  increase  the  reflectivity  of  the  specimen  so  that  they  could  better 
reflect  a  signal  back  to  the  laser  vibrometer,  the  beams  were  marked  with  a  white 
PENTEL  100W  S  paint  pen  in  the  manner  described  by  DeLeon  [1].  Ten  specimens 
were  tested  to  detennine  the  frequency  of  the  first  mode  of  the  beam  specimens.  Results 
of  the  test  are  shown  below  in  Table  3.1. 

Table  3.1:  Frequency  and  Mass  Experimental  Results. 


YSH-70-A/RS-3C  [0/90/0]  Specimens 


Sample 

Mode  1  (Hz) 

Mass  (mg) 

1 

150 

22.8 

2 

140.6 

23.5 

3 

140.24 

22.8 

4 

135.2 

22.7 

5 

148.4 

23.4 

6 

146.1 

23 

7 

150.8 

23 

8 

151.6 

23.7 

9 

137.5 

23.2 

10 

148.5 

24 

Average 

144.894 

23.21 

Standard  Variation 

5.9601911 

0.43063261 

A  FEA  model  was  constructed,  seen  in  Section  3.2.,  and  compared  against  the 
average  of  these  experimental  results. 


3.2.  FEA  Beam  Model 

An  FEA  model  of  the  carbon  fiber  test  specimens  was  constructed  using 
ABAQUS.  Two  models  were  made,  one  using  10  quadratic  (B32)  beam  elements  21 
nodes,  Figure  3.3,  and  the  other  using  64  composite  shell  elements  (SQR4),  264  nodes, 
Figure  3.4.  Both  models  used  the  material  properties  found  in  Section  2.2,  with  the 
necessary  calculations  and  changes  made  to  the  beam  element,  Chapter  2.  Note  the 
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I-beam  shape  of  the  beam  element  model  due  to  the  composite  structure.  The  base  of  the 
beams  were  clamped,  0  degrees  of  freedom  (DOF),  while  all  other  nodes  were  given  6 
DOF.  Modal  frequency  analyses  of  the  beam  elements  were  performed. 


Figure  3.3:  Beam  Element  Finite  Element  Model  of  Test  Specimens. 


Figure  3.4:  Composite  Shell  Element  model  of  Test  Specimens. 
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The  result  of  the  modal  analysis  of  both  beams  was  for  a  first  bending  frequency 
of  219.5  Hz,  Table  3.2.  This  was  significantly  lower  than  the  experimental  data  obtained 
in  Section  3.1.1.  Therefore,  the  composite  material  warranted  further  investigation  since 
the  same  material  and  manufacturing  process  was  being  used  in  the  engineered  wing.  It 
was  important  to  understand  what  was  causing  such  a  significant  deviation  in  the 
expected  modal  frequency  in  order  to  perfonn  an  accurate  analysis  of  the  wing. 


Table  3.2:  Comparison  of  Experimental  Data  to  FEA  Reasults. 


Comparison  of  Beam  Specimen  Results 

Experimental 
Beam  FEA 
Composite  Shell 

1st  Mode  (Hz)  %  Difference  to  FEA 

144.89  (avg)  34.0% 

219.5 

219.5 

3.2.1.  Modifiable  FEA  Beam  Model 

Since  there  was  a  significant  difference  between  the  experimental  and  analytical 
results,  a  MATLAB  code  was  developed,  Appendix  C.  This  code  imports  pre¬ 
determined  nodes  and  elements  based  on  the  input  file  (.inp)  of  the  wire  element 
developed  above  Appendix  D.  The  code  changes  the  thickness  of  the  specimen,  the  ply 
orientation  of  lamina,  and  the  angle  at  which  the  beam  was  cut  in  comparison  to  [0/90/0] 
orientation  of  the  composite  material.  These  changes  are  written  into  a  new  input  file  that 
is  solved  using  ABAQUS.  The  code  then  also  has  the  ability  to  read  the  data  file  (.dat) 
from  ABAQUS,  which  contains  the  modal  frequencies  of  the  model,  and  store  the  modal 
frequencies  in  MATLAB  where  they  can  easily  be  compared. 

This  code  allowed  rapid  changes  to  be  made  to  the  model,  therefore  allowing  a 
better  understanding  of  the  effects  of  the  thickness  of  the  specimen,  the  ply  orientation  of 
lamina,  and  the  angle  at  which  the  beam  was  cut  in  comparison  to  [0/90/0]  orientation  of 
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the  composite  material.  This  code  was  used  to  vary  each  of  these  inputs  to  allow  a  better 
understanding  of  how  they  affected  the  modal  frequency  of  the  test  specimens. 

3.3.  Carbon  Fiber  Analysis 

Several  options  were  explored  to  determine  the  source  of  error  in  the  laminate  specimens: 

1 .  Incorrect  material  properties  provided  by  the  manufacturer. 

2.  Incorrect  dimensions  due  to  low  number  of  plies  in  the  laminate,  and  small 
dimensions  of  the  specimens. 

3.  Errors  in  the  manufacturing  of  the  laminate  resulting  in  errors  in  ply  orientation. 
Each  of  these  was  examined  in  order  to  determine  a  root  cause  of  the  material  property 
deviation. 

The  first  thing  examined  were  the  given  material  properties  of  the  composite. 
O’Hara  had  previously  tested  twenty-ply  specimens  of  the  composite  using  ASTM  D 
3039  and  ASTM  D  3518  standards  shown  in  Table  2.1. 

3.3.1.  Laminate  Material  Dimensions 

Since  it  was  determined  that  the  material  properties  of  the  composite  matched 
closely  those  given  by  the  manufacturer,  the  next  thing  that  was  examined  was  the 
dimensions  of  the  test  specimens.  Incorrect  dimensions  could  cause  error  in  the 
equations  used  to  determine  the  expected  analytical  modal  frequency  of  the  test.  Using 
Fowler  Precision  calipers  with  a  tolerance  of  +/-  5  microns,  every  specimen  was  checked. 
Variation  observed  in  the  width  and  length  of  the  beam  observed  was  within  the  tolerance 
of  the  calipers  used,  and  therefore  was  not  investigated  further. 
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The  thicknesses  of  the  beams  were  also  measured.  It  was  found  that  the 


thickness  of  the  specimen  varied  between  135  and  160  microns,  and  had  an  average 
thickness  of  145  microns.  This  was  different  from  the  expected  150  microns,  and  has 
some  effect  in  the  modal  frequency.  A  finite  element  model  described  in  Section  3.2.1 
was  created  to  determine  how  much  of  an  effect  this  would  have  on  the  first  modal 


frequency.  Using  MATLAB  Code,  Appendix  C,  the  thickness  of  the  beam  was  varied 
from  135  to  160  microns,  resulting  in  6  runs  of  the  beam  model.  The  results  of  the  test 


are  shown  below  in  Figure  3.5. 


Modal  Freq.  V.  Thickness 


Beam  Specimen 


Figure  3.5:  Variation  of  Beam  Thickness  and  Effects  on  First  Modal  Frequency. 


The  results  of  the  test  shows  that  the  frequency  of  the  specimen  is  thickness 
dependent,  varying  between  197.56  and  234.09  Hz.  While  this  is  a  significant  change,  it 
does  not  explain  the  144.89  Hz  response  seen  during  the  experimentation.  Further 
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investigation  was  done  looking  at  the  ply  orientation  and  errors  that  could  have  been 
caused  during  manufacturing. 

3.3.2.  Lamina  Ply  Orientation  &  Manufacturing  Errors 

Errors  in  manufacturing  were  then  examined.  Based  on  the  calculation  in  Section 

2.2,  it  was  determined  that  a  small  variance  in  the  ply  orientation  could  cause  significant 
variations  on  the  material  properties.  The  alignment  of  the  fibers  was  then  examined  to 
see  if  such  a  variation  existed.  Multiple  test  specimens  measuring  40mm  long  by  1mm 
wide  were  constructed  out  of  the  same  [0/90/0]  carbon  fiber  composite  as  the  engineered 
wing  and  the  material  sample  specimens,  See  Figure  3.6.  From  these  specimens,  it  was 
hoped  to  determine  if  there  was  any  variation  seen  between  the  orientation  of  the  fibers 
and  along  the  straight  edge  of  the  test  specimen.  Any  such  variation  would  indicate  that 
either  the  ply  orientation  of  the  top  ply  was  not  at  a  0°  orientation,  the  angle  at  which  the 
laser  was  cutting  was  not  parallel  to  the  orientation  of  the  fibers,  or  both. 


Figure  3.6:  Test  Specimen  for  Ply-Orientation  Test. 

These  test  specimens  were  then  examined  using  a  Zeiss  Discovery  V.12  optical 
microscope,  (Figure  3.7).  This  microscope  has  up  to  150x  magnification,  allowing  the 
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individual  fibers  to  be  seen  on  the  composite  test  specimens.  Each  of  the  test  specimens 
were  examined  under  the  microscope,  with  digital  images  being  captured  for  each 
specimen,  (See  Figure  3.8).  The  fiber  orientations  of  the  test  specimens  were  then 
measured  using  a  developed  MATLAB  code,  (See  Appendix  E). 


Figure  3.7:  Zeiss  Discovery  V.12  Optical  Zoom  Microscope. 


Figure  3.8:  Image  of  Test  Specimen  used  to  Evaluate  Angle  of  Fibers  within  the 

Composite. 
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The  developed  MATLAB  code  is  designed  to  open  an  image  of  the  test  specimen 
in  MATLAB.  Once  in  MATLAB,  the  user  selects  two  points  on  the  straight  cut  edge  of 
the  test  specimen.  This  edge  was  selected  because  it  should  run  parallel  to  orientation  of 
the  fibers  on  both  the  top  and  bottom  plies  of  the  test  specimen.  The  code  then  forms  a 
line  between  the  two  points  to  use  as  a  reference  axis.  The  user  then  selects  two  points 
on  ten  separate  fibers  on  the  test  specimen.  Again,  the  code  forms  a  line  between  the  two 
points  selected  on  each  of  the  fibers.  After  the  points  on  all  ten  fibers  have  been  selected, 
the  code  calculates  the  average  angle  between  the  lines  formed  from  the  selected 
individual  fibers,  and  the  reference  axis.  This  process  was  repeated  several  times  for 
each  test  specimen  to  ensure  that  any  error  in  selecting  the  fibers  would  be  averaged  out. 


crrjtrd  with  the  HOI  li 


Selected  Points  (First  Click) 


Selected  Points  (Second  Click) 


Edge  (Reference)  Line 


Figure  3.9:  MATLAB  Code  Process  for  Measuring  Carbon  Fiber  Angle. 


58 


If  the  fibers  were  aligned  in  the  [0/90/0]  manner  as  prescribed  in  manufacturing, 
the  fibers  should  all  be  parallel  to  the  edge  of  the  test  specimen.  However,  this  test 
showed  that  this  was  not  the  case  for  any  of  the  specimens.  All  of  the  test  specimens 
exhibited  fibers  that  formed  a  small  angle  relative  to  the  reference  axis.  In  other  words, 
the  angle  of  the  fibers  were  not  at  [0/90/0]  degrees,  but  rather  had  some  degree  of 
variation  in  the  fiber  orientation.  Results  from  the  test  are  shown  below  in  Table  3.3. 

Table  3.3:  Results  from  Optical  Microscpe  Test. 


Optical  Microscope  Fiber  Measurements  1 

Test  Specimen 

Average  Angle 

1 

2.55° 

2 

3.49° 

3 

4.22° 

4 

3.08° 

5 

2.94° 

6 

2.33° 

7 

1.94° 

8 

1.65° 

9 

2.27° 

10 

2.75° 

Average: 

2.494° 

It  became  quickly  apparent  that  the  deviation  in  ply  orientation  was  the  result  of 
either  the  hand  lay-up  of  the  laminate,  or  the  placement  of  the  laminate  into  the  laser 
cutter.  Since  previous  testing  was  only  able  to  measure  the  angle  of  the  outer  two  layers 
relative  to  the  cut  of  the  specimen,  more  testing  was  done  to  confirm  the  sources  of  error. 

Wishing  to  quantify  if  a  variation  in  the  ply  angle  was  present,  two  other  sample 
carbon  fiber  specimen  tests  were  performed.  The  objective  of  the  tests  was  to  determine 
the  variation  in  the  orientation  of  the  carbon-fiber  lay-up  for  the  [0/90/0]  laminate,  as 
well  as  determine  the  ability  of  the  hand  lay-up  technique  to  correctly  achieve  the  desired 
ply  orientation. 
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Two  tests  were  devised  in  order  to  achieve  this  goal.  Both  test  specimens  were 
sheets  of  8”  x  11”  carbon  fiber.  The  composite  was  constructed  in  the  same  manner  the 
composite  used  in  the  previous  tests,  however  a  few  things  were  changed.  For  the  first 
test,  instead  of  all  three  plies  of  the  laminate  being  identical  sizes,  the  two  pre-preg  sheets 
fonned  the  middle  90°  layer  and  the  top  0°  layer  were  trimmed  to  be  slightly  smaller. 
The  pre-preg  used  for  the  top  layer  was  trimmed  slightly  smaller  than  the  pre-preg  used 
for  the  middle  ply.  As  the  laminas  were  stacked,  two  edges  of  the  lamina  were  aligned 
on  top  of  the  previous  ply.  This  meant  that  at  the  edge  of  the  laminate  sheet,  all  three 
plies  were  visible,  as  can  be  seen  in  Figure  3.10. 


Figure  3.10:  3-Ply  Test  Stack-up  Exposing  All  Three  Plies. 
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Figure  3.11:  3-Ply  Test  Specimen  Optical  Microscope  Digital  Image. 


Using  the  same  developed  MATLAB  code  as  the  previous  test,  the  angle  between 
the  plies  was  measured.  For  this  test,  the  bottom  ply  was  established  as  the  reference 
axis.  This  allowed  the  angle  of  the  other  fibers  to  be  measured  relative  to  that  bottom 
ply.  Figure  3.11  shows  an  image  taken  by  the  optical  microscope  of  the  test  specimen 
and  used  to  measure  the  angle  of  the  fibers. 

Another  test  was  performed  using  a  two-ply  0/0  laminate  to  again  verify  the 
results.  The  0/0  laminate  for  this  test  was  used  because  it  would  be  easier  to  visually  see 
the  error.  Again,  the  top  ply  was  trimmed  to  be  slightly  smaller  than  the  bottom  ply. 
Figure  3.12  shows  a  diagram  of  the  test  specimen.  Figure  3.12  shows  images  from  the 
optical  microscope.  The  same  develop  MATLAB  code  was  used  in  both  cases  to  evaluate 
the  image. 
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Figure  3.13:  Optical  Image  of  the  2-Ply  Unidirectional  Test  Specimen 
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The  results  from  both  tests  showed  that  the  same  approximate  2°  variations  in  ply 
angle  were  present  in  this  test  as  in  original  test.  Based  on  how  quickly  the  modulus  of 
the  lamina  properties  decrease  as  the  angle  of  the  ply  goes  away  from  0°,  the  error  in  the 
ply  orientation  was  primarily  responsible  for  the  low  first  modal  frequency  observed  in 
Section  3.1.1. 

3.4.  Sources  of  Error 

While  it  is  important  to  understand  the  effects  of  the  idiosyncrasies  in  the  ply 
orientation,  it  is  also  important  to  understand  how  they  became  present.  This  section  will 
serve  in  an  attempt  to  explain  why  these  occur. 

The  variation  in  thickness  of  the  ply  orientation  can  be  considered  an  effect  of 
nonnal  material  variation.  As  the  pre-preg  is  stacked,  there  is  some  bleed  off  of  the 
matrix  epoxy  as  it  is  cured  in  the  heat  press.  Since  the  laminate  contains  only  three  plies, 
even  a  small  amount  of  matrix  lost  correlates  to  a  larger  percentage  of  the  material  lost. 
While  a  5 -micron  loss  of  material  would  not  be  as  noticeable  in  a  laminate  with 
significantly  larger  number  of  plies,  the  thin  ply  nature  of  the  laminate  used  makes  it  a 
more  significant  factor. 

While  the  thickness  variation  can  be  considered  part  of  nonnal  material  variation, 
the  ply  orientation  cannot.  One  would  expect  that  the  ply  orientation  of  the  fibers  should 
vary  slightly  due  to  the  nature  of  the  hand  lay-up.  However,  the  fact  that  it  was 
consistently  off  needed  to  be  examined. 

One  of  the  things  assumed  about  the  carbon-fiber  pre-preg  sheets  was  that  they 
were  rectangular.  This  was  based  on  the  fact  that  the  roll  of  pre-preg  has  a  constant 
width,  and  all  cuts  were  made  parallel  to  the  edge  of  the  pre-preg  roll.  However,  it  was 
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found  that  the  edges  of  the  pre-preg  sheets  used  in  the  lay-up  of  the  composite  do  not 
always  remain  straight  during  construction,  Figure  3.14.  This  was  found  to  be  a  result  of 
the  material  moving  and  being  defonned  or  damaged  as  the  plastic  protective  sheet  was 
removed  on  the  pre-preg. 


Figure  3.14:  Image  showing  Defonnation  in  Pre-Preg  During  Manufacturing. 


Since  the  edges  of  the  composite  were  no  longer  straight,  the  use  of  a  straight 
edge  no  longer  ensured  that  the  angle  of  the  fibers  was  aligned  with  the  prescribed  fiber 
orientation,  Figure  3.15 
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Figure  3.15:  Small  Angle  Variation  in  Ply  Due  to  Deformed  Edge  of  Pre-Preg. 


Such  a  deformation  would  not  only  have  an  effect  as  the  pre-preg  plies  were 
stacked,  but  also  as  laminate  was  placed  into  the  laser  cutter.  Any  variation  in  the  angle 
in  which  the  composite  was  placed  into  the  laser  cutter  would  have  an  effect  on  the  ply 
orientation  by  essentially  “adding”  the  misaligned  angle  to  the  ply  orientation  of  each 
layer  of  the  composite,  Figure  3.16.  This  angle,  a,  either  can  be  caused  by  human  error, 
or  the  defonned  edge  of  the  composite  discussed  above. 


Figure  3.16:  Error  in  the  Placement  of  the  Composite  in  the  Laser  Cutter. 
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3.5.  Monte  Carlo  Validation 

This  following  section  shows  an  analytical  validation  of  the  experimental  results. 
These  results  are  important  to  show  the  range  at  which  the  variability  found  in  the 
experiments  can  have  on  the  effect  of  the  material  properties  and  modal  frequency  of 
carbon  fiber  specimens. 

A  Monte  Carlo  solution  to  determine  how  much  variation  in  ply  angle  (theta)  and 
laser  cutter  angle  (alpha),  and  thickness  would  account  for  the  differences  seen  in  the 
experimental  and  analytical  results.  A  Monte  Carlo  solution  uses  random  variables,  in 
this  case  using  a  Gaussian  distribution  of  random  variables,  to  vary  results.  The  Gaussian 
distribution  was  used  based  on  the  measurements  taken  of  the  various  variables.  The 
Monte  Carlo  approach  is  used  because  it  allows  a  quantification  of  the  variable  results 
seen,  as  well  as  to  show  the  probability  of  results  occurring. 

Using  the  same  FEA  model  developed  in  Section  3.2.1,  the  thickness  of  the 
specimen  (Mean  =  145  microns,  &  IStd),  the  ply  orientation  (Mean  =  2.5  deg,  &  1  Std) 
and  the  cut  angle  (Mean  =  2.5  deg,  &  1  Std)  were  varied  randomly  using  a  Gaussian 
distribution.  This  solution  was  done  using  ABAQUS  as  a  solver  and  developed 
MATLAB  code,  Appendix  F,  to  set  up  and  record  the  results  of  the  Monte  Carlo  solution. 
The  results  of  the  analysis  are  shown  in  Figure  3.17,  showing  the  first  modal  frequency 
for  each  of  the  runs  in  the  Monte  Carlo  solution. 
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Figure  3.17:  Results  of  Monte  Carlo  Solution. 
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Based  on  this  approach,  it  can  be  seen  that  from  the  observed  idiosyncrasies  in  the 
material,  the  low  first  modal  frequency  results  from  the  laser  testing  can  be  explained. 
Taking  into  account  all  the  idiosyncrasies,  an  average  first  modal  response  of  the 
specimens  tested  in  Section  3.1.1  is  154.16  Hz  according  to  the  Monte  Carlo  solution. 
Based  on  the  distribution,  the  144.89  Hz  average  can  be  explained  as  a  highly  probable, 
albeit  slightly  below  average,  solution. 

3.6.  Laminate  Material  Property  Analysis  Conclusions 

Based  on  the  experimentation  results  and  the  Monte  Carlo  analysis,  it  can  be  seen 

that  the  deviation  in  ply  orientation,  whether  caused  by  hand  lay-up  or  angle  in  which  the 
carbon  fiber  was  placed  into  the  laser  cutter,  and  the  slightly  smaller  thickness  were 
responsible  for  the  results  of  the  lower  first  mode  frequency  in  the  experimentation.  This 
could  lead  the  veins  of  the  engineered  wing  to  be  less  stiff  than  their  design  point,  and 
should  be  taken  into  account  when  perfonning  the  analysis  of  the  engineered  wing.  This 
also  underscores  the  importance  to  have  the  FEA  model  easily  changed  so  that  the 
variability  seen  with  the  composite  can  be  incorporated  when  comparing  the  results. 
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IV.  Manufactured  Wing  Experimentation  and  Analysis  Properties 

The  goal  of  this  project  was  to  analyze  the  manufactured  engineered  wing 
developed  by  O’Hara,  Figure  4.1.  The  veins  of  the  wing  were  constructed  according  to 
the  method  described  in  Section  2.1.1.  The  veins  of  the  wing  are  cut  out  of  a  single  sheet 
of  [0/90/0]  three-ply  unidirectional  carbon-fiber  composite,  based  on  the  geometry  of  the 
Manduca  Sexta,  with  a  Kapton  membrane. 


Vein  Structure  Membrane 


Figure  4. 1 :  Completed  Bio-Inspired  Wing. 

The  model  was  created  in  order  to  take  into  account  the  effects  of  the  carbon  fiber 
composite’s  idiosyncrasies,  allowing  for  the  idiosyncrasies  to  be  varied  within  the  model, 
similar  to  the  beam  analysis  perfonned  in  Section  3.5.  Both  analytical  and  experimental 
methods  were  used  in  order  to  evaluate  the  wing.  The  mode  shapes  and  frequencies  were 
used  to  compare  the  FEA  model  with  the  engineered  wing.  Norris  determined  that  the 
ratio  of  the  modal  response  of  the  wing  was  important  in  characterizing  the  wing’s 
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stiffness  [6].  Therefore,  it  was  important  the  FEA  model  matched  the  modal  frequency  of 
the  engineered  wing.  This  chapter  will  cover  the  set  up  of  the  experimentation,  the  FEA 
model  used  in  the  analysis.. 

4.1.  Engineered  Wing  Modal  Frequency 

DeLeon  measured  the  modal  shapes  and  frequencies  of  the  engineered  wing  [2]. 

It  is  these  experimental  values  for  which  the  model  was  set  up  and  compared  against.  For 
his  experiment,  the  base  of  the  wing  was  clamped,  Figure  4.2,  and  attached  to  a  Briiel  & 
Kjoer  Mini  Shaker  4810  that  would  impart  the  vibrations  into  the  wing. 


Figure  4.2:  Clamped  Wing. 


The  wing  was  placed  into  the  vacuum  chamber  at  AFIT,  Figure  4.3,  and  attached 
to  the  shaker.  A  pseudo  vacuum  (less  than  1%  atmospheric  pressure)  was  pulled,  and  the 
wing  was  shaken  using  a  pseudo  random  input.  The  dynamic  response  of  the  wing  was 
measured  using  a  Polytec  SLV.  The  SLV  uses  a  laser  to  take  precise  distance 
measurements,  and  is  capable  of  taking  measurements  along  several  points  along  the 
wing.  When  comparing  these  measurements  to  a  reference  point  on  the  clamp  of  the 
wing  and  the  pseudo  random  input  of  the  vibrometer,  the  SLV  is  able  to  determine  the 
modal  frequencies  and  corresponding  mode  shapes. 
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Figure  4.3:  Vacuum  Chamber  used  at  AFIT. 


Based  on  this  experiment,  DeLeon  determined  that  the  first  and  second  modal 
frequencies  of  the  engineered  wing  are  58.1  Hz  and  80.3  Hz  respectively.  These  values 
match  the  flap,  and  feather  mode  shapes  determined  by  Norris,  Figure  4.4.  [6]. 


Figure  4.4:  Mode  Shapes  of  the  Manduca  Sexta  Wing ,  flap,  feather,  saddle,  and  bisaddle. 
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4.2.  Wing  Model  MATLAB  Code 

One  of  the  main  purposes  of  this  project  was  to  create  an  FEA  model  of  the 
engineered  wing.  This  model  will  allow  for  any  changes  in  material  or  dimensional 
properties  of  the  wing  to  be  evaluated  prior  to  physical  fabrication,  saving  both  time  and 
money  in  the  creation  of  a  viable  FMAV.  Due  to  the  complex  nature  caused  by  the 
orientation  of  the  fibers,  and  the  curvature  of  the  veins,  Figure  4.5,  the  model  was 
developed  using  a  developed  MATLAB  code  to  generate  an  FEA  input  file  that  could  be 
solved  in  ABAQUS.  This  developed  code  made  it  easier  to  assign  individual  elements  on 
the  vein  correct  material  properties  based  on  the  curvature. 


Figure  4.5:  Diagram  Showing  Fiber  Orientation  over  the  Vein  Curvature. 

This  model  will  also  serve  to  evaluate  the  effects  of  the  carbon  fiber 
manufacturing  idiosyncrasies,  and  how  they  affect  the  modal  frequencies  of  the 
engineered  wing.  Therefore,  a  code  that  varied  the  idiosyncrasies  seen  within  the  carbon 
fiber  composite  was  set  up  and  an  analysis  was  performed  on  the  wing  model.  In  order 
for  this  approach  to  be  utilized,  the  model  needs  to  be  easily  modified  to  account  for  the 
varying  composite  material  properties.  Therefore,  a  MATLAB  code  was  developed, 
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Appendix  G,  so  that  multiple  input  files  can  be  created  to  reflect  such  varying 
manufactured  variations.  These  input  files  were  solved  using  ABAQUS,  with  the  results 
of  the  code  being  read  by  the  developed  MATLAB  code  to  collect  data  for  analysis.  This 
section  will  serve  to  explain  the  creation  and  capabilities  of  the  code  used  to  generate  the 
FEA  models.  The  flowchart  below,  Figure  4.6,  shows  the  general  process  perfonned  in 
the  creation  and  analysis  of  the  FEA  model. 

Furthermore,  since  the  experimental  modal  results  of  the  engineered  wing  did  not 
match  the  results  of  the  biological  wing,  the  model  was  made  to  be  easily  modified.  This 
would  allow  for  future  work  done  on  the  design  of  the  engineered  wing  to  be  easily 
incorporated  into  the  model.  The  model  has  the  ability  for  the  composite  ply  orientation, 
thickness  of  the  plies,  number  of  plies  and  width  of  the  veins  to  easily  be  changed.  Also, 
since  the  geometry  is  based  on  a  separate  file,  the  geometry  of  the  wing  can  easily  be 
changed  in  order  to  incorporate  camber  or  other  changes  of  geometry  compared  with 
using  ABAQUS’  CAD  program.  This  will  be  important  for  future  work  done  on  the 
wing. 
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Figure  4.6:  Flowchart  showing  the  Process  of  Generating  the  Wing  FEA  Model. 
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4.2.1.  Creating  the  Geometry  of  the  Wing 
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The  geometry  of  the  wing  was  created  in  MATLAB  using  a  code  developed  by 
O’Hara  [1],  The  code  bases  the  geometry  of  the  wing  on  an  image  taken  vertically  above 
a  Manduca  Sexta  wing  laid  on  a  flat  surface.  The  user  then  defines  the  points  along  the 
individual  veins,  which  are  then  splined  to  fonn  a  curve  that  matches  the  geometry  of  the 
vein,  Figure  4.7,  shown  in  blue. 


Spli 


Line  Formed  by 


of  Bio-Wing 


Figure  4.7:  Bio-wing  with  Selected  Points,  Creating  the  Basis  for  the  Vein  and 
Membrane  Nodes  for  the  Engineered  Wing  Model. 


Points  along  the  edge  of  the  membrane  are  then  selected,  Figure  4.7  (shown  in 
red).  The  program  then  creates  another  splined  line  along  the  outline  of  the  membrane. 
Using  the  points  created  along  the  lines  for  the  veins,  and  the  membrane,  various  sections 
of  membrane  were  created  that  could  easily  be  transfonned  by  ABAQUS  into  nodes  and 
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elements,  Figure  4.8.  The  generation  of  the  nodes  and  elements  will  be  explained  in 
Section  4.2.2.  The  geometry  was  visually  compared  to  the  engineered  wing  before 
proceeding  with  the  analysis  in  order  to  ensure  the  dimensions  matched. 


Membrane  Sections  Created 
by  Code  (Various  Colors) 


Figure  4.8:  Diagram  showing  membrane  sections  of  the  wing. 

4.2.2.  Generate  Nodes  &  Elements 

ABAQUS 

/ -  - V 

2.  Generates  Nodes  & 

Elements 

\ _ / 

The  lines  representing  the  veins  and  membrane  are  then  imported  into  ABAQUS, 
using  a  Python  script.  This  automatically  creates  points  with  X  and  Y  coordinates  along 
the  path  of  the  lines  created  with  O’Hara’s  program,  with  the  Z  coordinates  all  being  zero 
due  to  the  fact  that  the  wing  is  made  using  a  flat  plate.  These  points  will  serve  as  the 
location  for  the  nodes  within  the  model.  A  total  of  22,248  nodes  were  created  for  the 
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model.  Common  nodes  were  created  along  the  veins  which  were  used  for  both  the  veins 
and  the  membrane.  This  was  done  to  mimic  the  fact  that  the  membrane  was  physically 
attached  to  the  veins  using  adhesive. 

Elements  are  used  to  tie  the  nodes  of  the  model  together.  The  type  of  element 
used  affects  how  the  stiffness  matrix,  K,  of  the  model  is  set  up.  For  this  project, 
quadratic  beam  elements,  B32,  were  used  for  the  veins,  and  eight  node  shell,  S8R,  and 
six -node  shell,  STRI65,  elements  were  used  for  the  membrane. 

The  B32  element  is  a  3-node  quadratic,  1-D  beam  element  in  space,  Figure  4.9 
[22].  The  element  uses  parabolic  interpolation,  with  6  DOF,  3  translation,  about  the  X,Y 
and  Z  axis,  and  3  rotation,  about  the  X,Y,Z  axis.  This  Timoshenko  beam  element  was 
chosen  based  on  its  ability  to  predict  bending  displacements  compared  to  other  beam 
elements  [22]. 


The  S8R  shell  element  is  an  8-node  doubly  curved  shell  element  with  reduced 
integration,  Figure  4.10.  The  element  also  has  6  DOF  per  node  (3  translation,  3  rotation). 
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This  element  was  chosen  based  on  previous  experience  of  its  ability  to  model  the 
membrane  material.  The  S8R  element  does  not  have  the  ability  to  deviate  too  much 
from  its  square  shape  [23],  and  for  certain  cases,  which  will  be  explained  later,  triangular 
elements  were  used. 


Figure  4.10:  S8R  Element,  Showing  6  DOF. 


The  final  element  chosen  was  the  STRI65  element.  The  STRI65  element  is  a  6- 
node  element,  with  6  degrees  of  freedom  when  attached  to  6  DOF  nodes  on  certain 
elements  (such  as  S8R),  and  5  DOF  in  free  space,  or  at  boundary  conditions,  Figure  4.1 1 
[23].  Since  all  the  STRI65  elements  in  the  model  are  attached  to  S8R  elements,  6DOF 
exist.  Unlike  the  S8R  element,  this  element  does  not  have  the  ability  to  change  thickness, 
and  therefore  should  only  be  used  for  thin  membranes.  Since  the  membrane  of  the  wing 
is  thin,  20  microns,  relative  to  the  chordwise  direction,  50,000  microns  (50mm),  this  was 
an  acceptable  limitation.  This  element  was  also  chosen  based  on  its  ability  to  better 
model  certain  areas  within  the  geometry  due  to  angle  constraints. 
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Figure  4.1 1:  STRI65  Element. 


The  elements  were  assigned  to  the  nodes  using  a  build  in  ABAQUS  feature, 
where  the  user  would  select  the  geometry,  either  the  veins  or  the  membrane,  and  assign 
either  the  beam  or  the  plate  elements.  ABAQUS  used  both  triangular  and  quadrahedral 
elements  for  the  membrane  due  to  the  angle  restraints  placed  on  the  quadrahedral  shell 
elements.  Figure  4.12  shows  shell  elements  where  the  quadrilateral  formed  by  the  nodes 
has  interior  angles  that  are  within  acceptable  limits,  and  exceed  acceptable  limits  (greater 
than  145°,  less  than  45°)  where  a  triangular  element  would  be  used  instead. 


Nodes 


EA 


Good 


Bad 


Good 


Figure  4.12:  S8R  and  STRI65  Elements.. 
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The  vein  structure  of  the  wing  was  modeled  using  quadratic  beam  elements,  (B32 
elements  in  ABAQUS).  These  elements  represented  the  thin,  slender  shape  of  the  veins 
as  was  evident  from  the  material  property  evaluation  in  Chapter  III.  They  also  had  the 
advantage  over  the  use  of  shell  elements  in  that  they  could  be  easily  modified.  Unlike 
shell  elements,  the  width  of  veins  could  be  changed  easily  without  worrying  about  the 
geometric  considerations  of  the  four-node  elements,  such  as  the  angle  formed  between 
each  node  to  form  the  element,  Figure  4.13.  This  figure  shows  shell  elements  where  the 
quadrilateral  formed  by  the  nodes  has  interior  angles  that  are  within  acceptable  limits, 
and  exceed  acceptable  limits  (greater  than  135°,  or  less  than  45°). 


Nodes 


Good 


Figure  4.13:  Figure  Showing  both  Good  and  Bad  Placement  of  Nodes  for  Shell  Elements. 


If  the  geometry  of  the  wing  were  to  be  changed  for  future  iterations,  then  the  shell 
elements  would  have  to  be  checked  to  ensure  that  they  all  fall  within  acceptable 
limitations.  This  is  something  that  would  not  have  to  be  done  to  the  beam  elements  due 
to  the  fact  that  the  nodes  are  placed  in  a  straight  line  for  each  element.  The  width  of  the 
element  is  assigned  via  a  beam  section  profile.  The  profile  section  affects  the  stiffness 
matrix  of  the  model,  however  unlike  the  shell  element,  the  geometric  coordinates  of  the 
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nodes  and  elements  in  the  model  will  not  change,  Figure  4.14.  While  it  would  be 
possible  to  check  the  placement  of  the  nodes  in  the  shell  element  prior,  it  would  be  time 
intensive  and  provide  no  significant  advantage  versus  using  beam  elements  that  have 
already  proven  their  effectiveness. 


a)  Small  Profile  b)  Large  Profile 


Figure  4. 14:  Element  with  Two  Different  Profiles.  Notice  Actual  Element  and  Nodes 

Remain  Unchanged. 

1,093  B32  elements  were  created.  This  large  number  of  nodes  was  necessary  in 
order  to  ensure  that  angle  between  each  beam  element  was  small.  This  was  done  in  order 
to  better  capture  the  curved  geometry  of  the  wing,  as  well  as  the  material  properties, 
which  as  shown  previously  are  highly  sensitive  to  changes  in  ply  orientation.  1.093  is  the 
minimum  number  of  elements  that  will  keep  the  angle  between  each  beam  element  less 
than  one  degree,  except  for  where  the  Arculus  and  Cubitus  veins  meet  (see  Section 
4.2.4.). 

7,183  S8R  elements  and  175  STRI65  elements  were  used.  The  high  number  was 
necessary  to  ensure  that  the  membrane  could  be  fixed  to  the  veins  in  the  model,  i.e.  have 
common  nodes  for  the  beam  and  the  membrane.  Figure  4.15  shows  a  complete  image  of 
the  wing  showing  the  elements  used. 
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Figure  4.15:  Wing  Model  Showing  Elements. 


4.2.3.  Scale  Geometry 
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Since  the  geometry  in  O’Hara’s  program  is  outputted  in  inches,  the  scale  of  the 
geometry  needs  to  be  changed.  This  is  done  with  the  wing  model  MATLAB  code.  The 
code  determines  the  length  of  the  wing  based  on  the  X,Y,  and  Z  coordinates,  and  then 
scales  the  model  to  the  desired  size,  in  this  case  50mm.  This  is  important  because 
incorrect  geometry  can  have  a  significant  effect  on  the  modal  frequency  analysis. 

In  addition,  for  the  sake  of  simplicity,  the  model  was  shifted  so  that  every  node 
has  a  positive  X  coordinate.  This  was  done  in  order  to  simplify  calculations  done  later  in 
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the  program  dealing  with  the  beam  vein  element  orientation  angles  and  material 
properties. 

4.2.4.  Vein  Width 

MATLAB  Code 

( - \ 

4.  Sets  Vein  Widths 

V _ ) 

The  widths  of  the  individual  veins  of  the  manufactured  wing  were  varied  linearly 
from  root  to  tip  of  the  individual  vein.  This  was  done  as  an  approximation  based  on  the 
actual  Manduca  Sexta  wing  in  order  to  match  the  material  properties.  Figure  4.16  shows 
the  venation  patter  of  the  Manduca  Sexta  wing.  This  can  be  compared  to  the  geometry 
selected  for  the  engineered  wing,  Figure  4.17,  shown  in  blue.  In  order  to  simplify  the 
geometry  in  the  engineered  wing,  the  veins  Rl,  R2,  R3  and  R4,  and  A1  were  modeled  as 
one  vein. 
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Root 


Center  of  Pressure 


Costa  (C)  -  the  leading  edge  of  the  wing 

Subcosta(Sc)  -  second  longitudinal  vein  (behind  the  costa),  typically  unbranched 
Radius  (R)  -  third  longitudinal  vein,  one  five  branches  reach  the  wing  margin 
Media  (M)  -  fourth  longitudinal  vein,  one  to  four  branches  reach  the  wing  margin 
Cubitus  (Cu)  fifth  longitudinal  vein,  one  to  three  brances  reach  the  wing  margin 
Anal  Vein  (A)  -  unbranched  veins  behind  the  cubitus 

Figure  4.16:  Venation  Pattern  of  the  Manduca  Sexta. 


Figure  4. 17:  Venation  Pattern  of  the  Engineered  Wing. 


The  dimensions  of  the  veins  were  based  on  the  engineered  wing.  In  the  sizing  of 
the  engineered  wing,  the  width  of  the  veins  were  varied  linearly  from  the  root  of  the  vein 
to  the  tip.  For  manufacturing  simplicity,  the  Costal  and  Radial  veins  on  the  leading  edge 
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of  the  wing  were  combined  into  three  veins,  and  given  the  same  root  and  tip  dimensions, 
Figure  4.17.  Also,  the  Anal  vein,  the  Medial  veins,  and  two  of  the  Cubitus  veins  were 
grouped  (Figure  4.17)  and  given  the  same  dimensions  at  the  root  and  tip  of  the  veins. 
The  Cubitus  vein  along  the  Medial  Flexion  Line  (Figure  4.16  and  Figure  4.17),  and  the 
Arculus  were  given  unique  dimensions. 

Each  B32  element  in  the  above  stated  sets  of  veins  was  assigned  the  same  width 
for  the  root  and  tip.  Since  each  beam  element  must  maintain  a  constant  width,  the  value 
for  the  width  of  the  beam  element  was  detennined  based  on  the  location  of  the  mid-node. 
This  resulted  in  a  gradually  decreasing  beam  cross  section,  Figure  4. 18. 


4.2.5.  Node  and  Element  Sets 
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Node  sets  are  sets  of  nodes  that  either  are  assigned  a  boundary  condition,  or  are 
assigned  to  elements  that  are  part  of  the  same  element  set.  Element  sets  are  sets  of 
elements  that  share  that  same  material  property  and  section  property.  All  nodes  and 
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elements  in  the  model  must  be  assigned  to  a  node  set  and  an  element  set  in  order  to  be 
considered  in  the  stiffness  matrix  when  solving  the  model.  It  is  important  to  note  that 
this  is  not  where  any  properties  are  assigned  to  the  nodes  or  elements,  rather  the  nodes 
and  elements  are  just  being  sorted  into  groups  to  be  assigned  those  properties  later  in  the 
code. 

Assigning  the  nodes  and  elements  to  the  various  node  and  element  sets  was  done 
in  the  wing  model  MATLAB  code.  In  the  case  of  this  model,  the  nodes  that  make  up 
each  beam  element  need  to  be  placed  in  their  individual  node  set.  This  is  because  as  the 
curvature  of  the  beam  changes,  so  do  the  material  properties,  and  beam  profile,  See 
Section  4.2.7.  This  resulted  in  1,093  different  node  sets  being  created  to  account  for  the 
beam  elements.  Three  more  sets  were  created,  one  for  the  S8R  elements,  one  for  the 
STRI65  elements,  and  another  for  the  boundary  condition.  This  resulted  in  a  total  of 
1,096  node  sets.  Only  1,095  element  sets  were  created  due  to  the  fact  that  the  boundary 
condition  is  only  applied  to  nodes,  and  not  element.  Therefore,  there  was  no  need  to 
create  an  element  set  for  the  boundary  conditions. 
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4.2.6.  Determine  Material  and  Profile  Properties 


Repeat  for  Each 
Beam  Section 


MATLAB  Code 


Material  Property 


The  material  properties,  and  profile  section  of  the  elements  of  the  model  play  an 
important  part  in  determining  the  stiffness  of  the  model.  These  properties  need  to  be 
defined  for  every  element  in  the  model  in  order  to  form  the  stiffness  matrix,  K.  Since  the 
Kapton  membrane  is  an  isotropic  material,  all  the  shell  elements  representing  the 
membrane  can  share  the  same  material  property  and  profile  section.  However,  due  to  the 
changing  local  axis  of  the  material  due  to  the  curvature  of  the  veins,  Figure  4.19,  each 
B32  element  required  its  own  unique  material  property  and  beam  section  based  on  the  ply 
orientation  at  the  local  axis.  This  section  will  discuss  the  process  used  to  determine  and 
assign  the  section  properties  and  material  properties  for  each  set  of  elements  created  in 
the  previous  section. 
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Figure  4.19:  The  Engineered  Wing  Showing  the  Directions  of  the  Composite  Fibers. 


Kapton’s  properties,  Table  4.1,  were  applied  to  all  shell  S8R  and  STRI65 
elements  using  a  single  set  of  material  properties. 


Table  4.1:  Kapton's  Material  Properties 


Material  Properties,  Kapton 

Modulus 

2.5  G  Pa 

Density 

1.42  g/cc 

Poisson’s  Ratio 

0.3 

Thickness 

20  microns 

The  application  of  the  beam,  B32  elements  properties  required  more  effort  due  to 
the  anisotropic  nature  of  the  material.  Since  the  material  property  varied  based  on  the  ply 
angle,  the  local  axis  needed  to  be  determined.  Therefore,  based  on  the  coordinate 
location  of  the  two  end  nodes  in  each  element,  the  angle,  which  the  element  made  in 
relation  to  the  global  axis,  was  determined  in  the  MATLAB  code  using  trigonometry. 
This  was  done  in  order  to  create  a  homogeneous  beam,  which  could  more  easily  be 
modified  than  the  ABAQUS  composite  shell  element,  especially  when  changing  the 
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physical  width  of  the  veins.  This  avoided  problems  with  shell  element  nodal  angles  as 
described  in  Section  4.2.2. 


Figure  4.20:  Difference  betweeen  the  Global  and  Local  Axis  for  the  Beam  Element. 

Once  the  angle  for  each  element  has  been  detennined,  the  local  axis  for  the 
material  can  be  determined.  This  is  done  by  subtracting  the  difference  between  the  angle 
fonned  between  the  local  and  the  global  axis.  This  was  done  for  each  ply  angle  in  the 
[0/90/0]  composite  and  resulted  in  the  ply  orientation  angle  for  each  lamina.  This 
infonnation  was  necessary  in  order  to  determine  the  material  properties  of  the  composite 
for  the  element. 

After  the  local  ply  orientation  for  the  element  was  determine,  the  local  material 
properties  were  determined,  based  on  Equation  2.1  in  Section  2.3,  using  a  material 
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property  MATLAB  function,  Appendix  B.  This  code  determined  the  local  lamina 
properties,  which  would  be  used  during  the  generation  of  the  beam  profde  cross  section. 

In  order  to  account  for  the  differences  in  the  modulus  for  the  individual  lamina, 
the  effective  moment  of  inertia  of  the  element  needed  to  be  detennined.  The  MATLAB 
code  then  perfonned  the  effective  moment  of  inertia  calculations  described  in  Section 
2.5.2.,  Figure  4.21.  The  width  of  the  vein  was  detennined  based  on  the  input  given  in 
Section  4.2.4.  In  order  to  keep  the  vein  visualization  in  ABAQUS  consistent  with  the 
physical  widths  of  the  veins  at  each  element  location,  the  largest  modulus,  Ex,  calculated 
amongst  the  three  plies  was  used  to  determine  the  effective  width  of  the  other  two  plies. 


Figure  4.21 :  Effective  Moment  of  Inertia  Transformation  for  the  Vein  B32  Profile. 

This  method  of  transformation  was  only  applicable  as  long  as  effective  width  of 
the  top  on  bottom  lamina  was  greater  than  Vi  the  effective  width  of  the  middle  lamina, 
See  Figure  4.22.  ABAQUS  was  unable  to  perform  the  calculations  on  profile  sections 
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that  exceeded  this  criterion  because  it  is  incapable  of  solving  means  where  the  moment  of 
inertia  in  the  2  direction  is  larger  than  the  moment  of  inertia  in  the  1  direction. 


Good 


Good 


Good 


Bad 


Figure  4.  22:  Effective  Moment  of  Inertia  Profile  Sections. 


Six  elements  in  the  original  model  failed  this  criterion,  all  of  which  were  in  the 
Arculus  where  the  angle  of  the  beam  element  approaches  90°  from  the  global 
coordinates.  In  order  to  deal  with  these  six  elements,  a  check  was  put  in  the  code  that 
would  transfonn  bad  elements  into  rectangular  beams  with  the  same  profile  cross 
sectional  area,  Figure  4.23.  Since  this  was  only  done  in  for  six  of  the  1,093  beam 
elements,  the  effect  in  the  modal  frequency  was  assumed  to  be  negligible. 


_ n _  ! - 

u  - 

Bad  Good 

Figure  4.23:  Changing  of  Cross  Sectional  Profile  for  Beam  Profiles  at  Angles  Close  to 

90°. 


As  the  profiles  were  formed,  the  material  properties  of  the  largest  modulus,  were 
recorded  and  placed  into  the  input  file  as  material  properties.  A  total  of  1,093  material 
properties  were  calculated,  one  for  each  beam  element.  These  properties  were  recorded 
with  the  material  property  for  the  Kapton  in  the  generated  input  file.  Figure  4.24  shows 
an  up-close  view  of  a  vein  and  the  vein  profile  used  in  the  wing. 
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Figure  4.24:  Wing  Model  Elements  Zoom-in  View,  Showing  Elements  and  Beam  Profile. 


4.2. 7.  Boundary  Conditions 

MATLAB  Code 

/ -  - \ 

7.  Specifies  Model's 

Boundary  Conditions 

V _ ) 

A  boundary  condition  in  finite  element  method  is  used  to  limit  the  degrees  of 
freedom  of  each  node.  In  the  case  of  the  project,  the  boundary  condition  present  was  the 
clamped  base  of  the  wing,  Figure  4.25.  For  this,  a  node  set  was  creating  that  included  6 
nodes  at  the  base  of  the  model  that  were  limited  to  0  degrees  of  freedom,  i.e.  no 
translation  or  rotation.  Figure  4.26  shows  the  location  of  the  clamped  nodes  on  the  FEA 
model. 
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Figure  4.25:  Clamped  Base  of  the  wing. 


Figure  4.26:  Wing  Model  Showing  Clamped  Area. 

4.2.8.  Output  Requests 


MATLAB  Code 


8.  Specifies  ABAQUS 
OUTPUT  (MODAL 
FREQ) 

V _ > 


A  step  module  in  ABAQUS  define  the  type  of  analysis  that  will  be  performed  on 
the  model.  A  step  also  defines  any  loads  and  boundary  conditions  that  will  be  performed 
in  the  analysis.  For  the  case  of  this  project,  the  clamped  boundary  condition  described  in 
Section  4.2.7.  was  applied.  A  step  was  created  in  order  to  solve  the  modal  frequency  of 
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the  model.  The  first  10  modes  of  the  model  were  solved  for  in  the  solution.  No  loads 


were  applied  due  to  the  nature  of  modal  frequency  analysis. 


4.2. 9.  Processing  Data 


Multiple  Run  Code 


11.  Multiple  Runs 
Performed 


ABAQUS 


~T 


i 


Creates  Text  File 
_ _ _ ^ 


Generated  Input  File 


Data  File 


MATLAB  Code 

^  N 

4.  Sets  Vein  Widths 

I 


10.  Processes  Data 
From  ABAQUS 


Modal  Frequency  Data 


Once  this  step  is  complete,  the  creation  of  the  model  is  complete,  and  it  is  ready  to 
be  solved.  The  MATLAB  code  develops  an  input  file,  Appendix  H,  based  on  the 
conditions  specified  in  this  chapter.  Figure  4.27  shows  a  completed  model  of  the  wing 
with  the  elements  shown.  The  beam  profiles  are  displayed  and  shown  in  green,  with  the 
membrane  being  displayed  in  white.  The  input  file  needed  to  be  solved  using  a  finite 
element  solver,  and  multiple  runs  needed  to  be  conducted  to  determine  the  effect  of  the 
carbon  fiber  idiosyncrasies  on  the  modal  frequencies  of  the  wing  will  be  discussed. 
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Figure  4.27:  Final  Model  of  the  Wing  with  Beam  Elements  Highlighted(Green). 


ABAQUS  was  used  in  order  to  solve  the  developed  model.  This  was  done 
automatically  using  MATLAB  commands,  without  manually  importing  the  file  into 
ABAQUS  CAE  itself: 

(eval ( [ ' dos ( ' ' abq6111  job='  Input_File  '  interactive'')']))-  (Note,  the 
.inp  file  must  be  stored  in  the  open  MATLAB  directory). 

ABAQUS  creates  a  data  file  with  the  outputs,  eigenvalues  and  eigenvectors, 
specified  in  Section  4.2.8.  The  data  file  is  a  text  file,  and  is  read  by  the  MATLAB  code, 
which  stores  the  values  for  the  first  two  modes,  which  were  visually  determined  to  be  the 
flap  and  feather  modes,  in  a  matrix  to  be  graphed  and  analyzed. 

The  Multiple  Run  Code  code  was  developed,  Appendix  I,  in  order  to  submit 
variable  runs  to  ABAUS  in  order  to  determine  the  effects  of  the  carbon  fiber 
idiosyncrasies  on  the  wing.  As  with  the  beam,  it  varied  the  thickness,  the  ply  orientation 
and  the  cut  angle  of  the  vein  structure  of  the  wing.  However,  due  to  the  time  constraints 
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and  the  significantly  longer  run  time  for  the  wing  analysis  as  opposed  to  the  beam,  a  full 
Monte  Carlo  solution  was  deemed  unrealistic,  as  it  would  have  taken  months  to  complete 
the  number  of  runs  required  for  a  full  solution.  Instead,  the  effects  of  the  individual 
variables,  thickness,  laser  cut  angle,  and  ply  orientation  were  varied  to  detennine  the 
probable  range  of  modal  frequencies  that  the  wing  could  experience. 
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V.  Manufactured  Wing  Analysis  Results  and  Discussion 


The  goal  of  this  chapter  is  to  discuss  and  evaluate  the  results  of  the 
experimentations  performed  on  the  engineered  wing  in  Chapter  4.  The  results  will  first 
cover  an  ideal  solution  that  does  not  take  into  account  any  variations  in  the  composite 
material  discussed  in  Chapter  3,  and  second,  multiple  runs  of  the  FEA  analysis  were 
perfonned  where  the  effects  of  the  variable  factors,  thickness,  theta  orientation,  9,  and 
placement  in  the  laser  cutter,  a,  were  tested  independently  in  order  to  compare  their 
effects  on  the  first  and  second  modal  frequencies  of  the  model.  Finally  using  the 
information,  the  experimental  modal  frequencies  determined  were  matched  based  on 
observed  idiosyncrasies  in  the  composite  material. 

5.1.  Ideal  Engineered  Wing 

The  first  model  solved  was  of  the  ideal  engineered  wing.  This  FEA  model  did  not 
take  any  of  the  variations  observed  in  Chapter  3  into  account,  and  represents  the 
engineered  wing  as  designed.  The  carbon  composite  in  the  model  was  set  to  a  ply 
orientation  of  [0  /90/0],  with  a  cut  angle,  a,  of  0°,  and  a  thickness  of  150  microns.  The 
calculated  mass  of  the  model  is  shown  in  Table  5.1.  The  modal  frequencies  of  the  first 
four  mode  shapes  are  listed  in  Table  5.2. 

Table  5.1:  Calculated  Mass  of  the  Ideal  Engineered  Wing. 


Mass  of  the  Wing 


FEA  Experimental  Difference 
Mass  51.3  mg  52.5  mg  2.3% 
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Table  5.2:  Modal  Frequencies  for  the  Ideal  Engineered  Wing. 


Modal  Frequencies  of  the  Ideal  Engineered  Wing 

Mode 

FEA  Frea.  [Hzl 

Experimental  Frea.  [HzJ 

Difference 

1 

62.6 

58.1 

7.19% 

2 

73.0 

80.3 

10.0% 

3 

138.6 

- 

- 

4 

198.5 

- 

- 

The  first  four  mode  shapes  are  shown  in  Figures  5.1  through  5.4.  Additional 
images  are  available  in  Appendix  J. 


Figure  5.1:  First  Mode  Shape  of  the  Engineered  Wing. 
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Figure  5.2:  Second  Mode  Shape  of  the  Engineerd  Wing. 


Figure  5.3:  Third  Mode  Shape  of  the  Engineered  Wing. 
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Figure  5.4:  Fourth  Mode  Shape  of  the  Engineered  Wing. 

The  mode  shapes  determined  by  the  FEA  model  closely  match  the  flap,  feather, 
saddle  and  bisaddle  shapes  exhibited  by  the  biological  Manduca  Sexta  wing  and  those 
determined  experimentally  of  the  engineered  wing  by  DeLeon,  Figure  5.5  [2].  It  can 
also  be  seen  that  the  mass  of  the  model  calculated  by  ABAQUS  is  within  2.3%  of  the 
experimentally  calculated  mass.  The  difference  of  1.2  mg  can  be  accounted  for  due  to 
tolerance  of  the  scale,  or  in  the  application  of  the  adhesive  for  the  Kapton  membrane. 


Figure  5.5:  First  Four  Mode  Shapes  of  the  Hawkmoth  Wing. 
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However,  there  is  a  considerable  difference  of  7.19%  and  10.0%  between  the  first 
and  second  modes  respectively  of  the  FEA  results  with  the  experimental  results.  Based 
on  the  experimental  results  of  the  composite  beam  structure,  Chapter  3,  it  was  considered 
highly  probable  that  the  variations  present  in  the  composite  material  could  account  for  the 
difference  between  the  experimental  and  FEA  results.  Therefore,  the  effects  of  the 
carbon  fiber  composite  variations  (idiosyncrasies)  were  examined  further. 

5.2.  Effect  of  Individual  Variables 

In  order  to  understand  the  effects  the  thickness,  cut  angle,  a,  and  ply  orientation, 
9,  have  on  the  dynamic  response  of  the  wing,  each  of  the  variables  was  varied 
independently  while  the  other  variables  were  held  at  the  initial  ideal  conditions.  The  first 
and  second  modal  frequencies  were  then  solved  for  and  compared  to  the  ideal  case  in 
order  to  detennine  the  effects  of  each  of  the  variables  on  the  wing.  Unlike  with  the  beam, 
the  angles  were  varied  in  both  the  positive  and  negative  direction  since  a  positive  and 
negative  angle  of  the  same  magnitude  could  have  different  effects  that  would  not  be 
present  on  the  symmetric  beam. 

5.2.1.  Thickness  of  the  Composite  Vein  Structure 

The  thickness  of  the  composite  veins  was  varied  in  the  FEA  model  between  the 
observed  135  micron  to  165  micron  variation.  The  cut  angle  and  ply  orientation  were 
held  constant  at  0,  and  [0/90/0]  respectively.  The  results  are  shown  in  Figure  5.6,  and 
Table  5.3  below.  For  each  run,  the  first  and  second  mode  shape  was  visually  determined 
to  match  the  flap  and  feather  modes  exhibited  in  the  ideal  case. 
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Table  5.3:  Effect  of  Thickness  on  Modal  Frequency. 


Effect  of  Thickness  on  Modal  Frequencies  of  the  Engineered  Wing 

Thickness  ( microns ) 

First  Modal  Frea.  (Hz) 

Second  Modal  Frea .  (Hz) 

135 

55.6 

66.1 

140 

58.0 

68.4 

145 

60.3 

70.7 

150 

62.6 

73.0 

155 

65.0 

75.4 

160 

67.3 

77.7 

165 

69.6 

80.0 

The  results  show  that  the  thickness  of  the  composite  has  a  significant  role  in 
determining  the  modal  frequency  of  the  wing.  While  the  mode  shapes  retain  the  same 
shape,  the  first  modal  frequency  varies  between  55.68  Hz  and  69.59  Hz,  while  the  second 
modal  frequency  varies  between  66. 13  Hz  and  79.99  Hz.  Both  modes  increase  linearly  as 
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the  thickness  of  the  composite  is  increased.  The  equations  and  R-values  are  shown  in 
Figure  5.6.  The  upward  trend  in  frequency  is  expected,  as  a  thicker  beam  member  would 
have  a  higher  moment  of  inertia,  and  therefore  be  more  resistant  to  bending. 


5.2.2.  Ply  Orientation  Angles,  Theta 

For  the  variation  in  the  ply  angle,  the  angle  of  the  mid-ply,  and  the  angle  of  the 
top  and  bottom  plies  were  varied  in  two  different  set  of  FEA  runs.  This  was  done  to 
show  the  effect  of  the  0°  plies  and  the  90°  plies  independently.  Due  to  quasi  symmetry  of 
the  actual  ply,  both  the  top  and  bottom  plies  were  varied  together  during  the  analysis 

The  mid-ply  of  the  FEA  model  was  varied  from  85°  to  95°.  The  results  of  the 
FEA  are  shown  below  in  Figure  5.7  and  Table  5.4.  The  thickness  of  the  composite  was 
held  constant  at  150  microns,  and  the  cut  angle  was  set  to  0°  for  each  of  the  FEA  runs. 


Effect  of  Middle  Ply  Orientation  on  Modal 

Freq. 


-S  63 
o 

^  61 


59 

57 

55 

84  86  88  90  92  94  96 

Alpha  (deg.) 


First  Mode 
•  Second  Mode 


Figure  5.7:  Effect  of  the  Mid-Ply  Orientation  of  the  Wing. 
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Table  5.4:  Effect  of  Mid-Ply  Angle  on  Modal  Frequencies. 


Effect  of  Mid-Ply  Orientation  on  Wing  Modal  Frequencies 

Plv  Orientation  (°) 

First  Modal  Freq.  (Hz) 

Second  Modal  Freq.  (Hz) 

[0/85/0] 

62.6 

72.9 

[0/86/0] 

62.6 

72.9 

[0/87/0] 

62.6 

73.0 

[0/88/0] 

62.6 

73.0 

[0/89/0] 

62.6 

73.0 

[0/90/0] 

62.6 

73.0 

[0/91/0] 

62.6 

73.1 

[0/92/0] 

62.7 

73.1 

[0/93/0] 

62.7 

73.2 

[0/94/0] 

62.7 

73.2 

[0/95/0] 

62.7 

73.2 

As  can  be  seen  in  Figure  5.7,  the  orientation  of  the  mid-ply  does  not  have  a 
significant  effect  on  the  modal  frequency  of  the  wing.  The  first  modal  frequency 
increases  from  62.647  Hz  to  62.655  Hz.  The  second  mode  increases  from  72.887  Hz  to 
73.244  Hz.  This  is  a  maximum  of  a  0.00013%  increase  in  the  first  modal  frequency,  and 
maximum  of  a  0.49%  increase  in  the  second  modal  Frequency.  Variances  in  modal 
frequency  this  small  would  be  hard  to  quantify  experimentally,  and  does  not  contribute 
significantly  to  the  variation  measured  by  DeFeon  and  the  ideal  wing. 

Another  test  was  performed  on  the  top  and  bottom  plies  of  the  composite  vein 
structure.  The  ply  angle  for  these  lamina  was  varied  from  -5°  to  5°  in  the  FEA  model. 
The  laser  cut  angle  was  held  constant  at  0°,  and  the  thickness  constant  at  150  microns. 
Figure  5.8  shows  the  results  of  the  analysis. 
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Effect  of  Top  and  Bottom  Ply  Orientation  on 

Modal  Freq. 


-6  -4  -2  0  2  4  6 

Theta,  Top  and  Bottom  ply  (deg.) 


First  Mode 
•Second  Mode 


Figure  5.8:  Effect  of  Top  and  Bottom  Ply  Orientation  on  the  Wing. 


Table  5.5:  Effect  of  Top  and  Bottom  Ply  Orientation  on  Wing  Modal  Frequencies 


Effect  of  Top  and  Bottom  Ply  Orientation  on  Wing  Modal  Frequencies 

Plv  Orientation  (°) 

First  Modal  Freq.  (Hz) 

Second  Modal  Frea.  (Hz) 

[-5/90/-5] 

53.7 

77.2 

[-4/90/-4] 

55.5 

76.6 

[-3/90/-3] 

57.3 

75.9 

[-2/90/-2] 

59.2 

75.0 

[-1/90/-1] 

61.0 

74.0 

[0/90/0] 

62.6 

73.0 

[1/90/1] 

64.0 

72.2 

[2/90/2] 

64.6 

71.8 

[3/90/3] 

64.5 

72.0 

[4/90/4] 

63.6 

72.6 

[5/90/5] 

62.3 

73.3 

Unlike  the  orientation  of  the  mid-ply,  the  orientation  of  the  top  and  bottom  ply 
has  a  much  more  significant  effect  on  the  modal  frequency  of  the  wing.  The  first  modal 
frequency  varies  from  53.719  Hz  at  a  9  of  [-5/90/-5]  to  64.633  Hz  at  a  9  of  [2/99/2].  The 
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results  form  a  concave  down  curve,  peaking  at  approximately  2°.  The  opposite  effect  of 
the  ply  orientation  can  be  seen  on  the  second  modal  frequency.  The  second  modal 
frequency  varies  from  a  high  of  77.199  Hz  at  a  9  of  [-5/90/-5]  to  a  low  of  71.846  Hz  at  a 
9  of  [2/992].  The  results  form  a  concave  up  curve  with  a  low  point  of  approximately  2°. 
The  difference  of  19.91  Hz  and  5.353  Hz,  seen  on  the  first  and  second  modal  frequencies 
respectively,  needed  to  be  examined  further. 

These  results  also  show  that  the  ply  orientation  of  the  top  and  bottom  plies  not 
only  play  an  important  part  in  determining  the  modal  frequency  of  the  wing,  but  unlike 
the  variation  in  the  thickness  of  the  composite,  it  also  plays  an  important  role  in 
detennining  the  modal  ratio  (MR,  the  ratio  between  the  first  and  second  modal 
frequencies)  of  the  wing.  As  the  orientation  of  the  ply  deviate  from  approximately  2°,  the 
MR  of  the  wing  increased.  This  is  due  to  the  fact  that  as  the  orientation  of  the  fibers 
deviate  from  this  2°  orientation  ,  the  equations  in  Chapter  2  show  that  the  fibers  will 
have  a  higher  elastic  modulus  in  the  chordwise  direction,  and  a  lower  modulus  in  the 
spanwise  direction.  This  would  make  the  wing  less  stiff  for  the  flap  mode,  thereby 
reducing  the  modal  frequency  of  that  mode,  and  stiffer  against  the  feather  mode,  a  torsion 
mode  in  the  chordwise  direction. 

The  curvature  of  the  veins  has  a  significant  effect  on  the  modal  frequency.  Since 
all  of  FEA  models  were  run  with  the  same  geometry,  the  fiber  angles  were  examined  and 
compared  to  the  geometry  of  the  wing,  Figure  5.9. 
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Figure  5.9:  Vein  Structure  of  the  Wing  Compared  to  Positive  and  Negative  Ply 

Orientation. 

As  Figure  5.9  shows,  as  the  ply  orientation  deviates  from  0°  in  the  negative 
direction,  the  fibers  tend  to  run  closer  to  parallel  along  the  Cubitus  vein,  as  they  deviate 
in  the  positive  direction,  they  tend  to  run  closer  to  parallel  along  the  leading  edge  veins. 
The  closer  the  fibers  are  to  running  0°  to  the  local  axis  along  either  of  these  veins,  the 
stiffer  these  veins  become.  Since  it  has  been  demonstrated  how  significantly  the  material 
properties  of  the  veins  change  as  the  fibers  approach  0°  to  the  local  axis,  any  positive  or 
negative  ply  orientation  would  have  a  significant  effect  on  the  stiffness  of  these  veins. 

Based  on  these  results,  it  can  be  concluded  that  the  stiffness  of  these  veins  play  an 
important  part  in  determining  the  modal  frequencies  of  the  wing.  In  addition,  since  the 
fibers  are  unidirectional,  as  the  fibers  are  orientated  to  stiffen  one  of  these  veins,  the 
effect  would  inversely  affect  the  other. 
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5.2.3.  Laser  Cut  Angle,  Alpha 


The  laser  cut  angle,  a,  of  the  composite  (the  angle  in  which  the  composite  in 
placed  in  the  laser  cutter  compared  with  the  desired  orientation,  Figure  5.10)  was  varied 
in  the  FEA  model  between  -5°  and  5°.  The  thickness  and  ply  orientation  were  held 
constant  at  150  microns,  and  [0/90/0]  respectively.  The  results  displayed  in  Figure  5.11 
and  Table  5.6  below. 


108 


Effect  of  Cut  Angle  on  Modal  Freq. 
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Alpha  (deg.) 


First  Mode 
Second  Mode 


Figure  5.11:  Effect  of  the  Laser  Cut  Angle  on  the  Wing 


Table  5.6:  Effect  Of  Laser  Cut  Angle  on  Engineered  Wing  Modal  Frequencies. 


Effect  of  Laser  Cut  Angle  on  Wing  Modal  Frequencies 


Laser  Cut  Angle  (°) 

-5 

-4 

-3 

-2 

-1 

0 


First  Modal  Freq.  (Hz) 

53.7 

55.5 

57.3 
59.2 
61.0 

62.6 
64.0 

64.7 
64.5 

63.7 

62.4 


Second  Modal  Freq.  (Hz) 

77.2 

76.6 

75.8 
75.0 
74.0 
73.0 

72.3 

71.9 
72.0 

72.6 
73.2 


As  can  be  seen  in  Figure  5.11,  the  effect  of  the  cut  angle  has  a  similar  effect  as  the 
effect  of  the  top  and  bottom  plies’  orientation.  There  is  less  than  a  0.1  Hz  difference 
between  the  results  of  the  test  that  varied  the  top  and  bottom  plies’  orientation  Table  5.5, 
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and  the  test  that  varied  the  laser  cut  angle.  The  laser  cut  angle  simply  adds  or  subtracts 
the  angle  to  the  ply  orientation.  Since  the  variation  in  the  orientation  of  the  mid  ply  has 
very  little  effect  on  the  modal  response  of  the  wing,  it  can  be  concluded  that  the  effect  of 
the  laser  cut  angle  is  simply  having  the  same  effect  as  the  fiber  orientations  of  the  top  and 
bottom  plies. 

5.3.  Effect  of  Composite  Thickness  and  Fiber  Angle 

Since  the  thickness  and  fiber  angle  of  the  top  and  bottom  fibers  had  the  most 

effect  on  the  frequency  response  of  the  wing,  FEA  runs  were  perfonned  that  varied  both 
the  thickness  and  laser  cut  angle  of  the  composite  for  the  wing.  The  laser  cut  angle  was 
chosen  over  varying  the  top  and  bottom  plies’  orientation  because  it  has  the  same  effect 
as  varying  the  top  and  bottom  ply  orientations,  would  be  simpler  to  modify  within  the 
program,  and  because  the  effect  of  the  variation  in  the  mid-ply  was  deemed  negligible. 

In  order  to  detennine  the  range  of  frequencies  that  could  be  expected  within  the 
wing,  the  laser  cut  angle  was  varied  between  -6°  and  6°  in  increments  of  1°.  Since  the 
thickness  caused  both  the  first  and  second  modal  frequency  to  vary  linearly,  the  values 
for  the  thickness  were  only  varied  using  the  minimum  and  maximum  values  of  135  and 
165  microns.  Again,  the  results  were  visually  checked  to  ensure  that  the  mode  shapes 
matched  those  exhibited  in  the  ideal  case,  Section  5.1.  The  results  for  the  first  modal 
frequency  is  shown  in  Figure  5.12,  and  the  second  modal  frequency  in  Figure  5.13,  both 
results  are  shown  in  Table  5.7 


110 


Modal  Freq.  vs.  Overall  Ply  Orientation  Angle 


-8  -6  -4  -2  0  2  4  6  8 


—♦—Mode  1, 135 
Microns 


Mode  1, 
microns 


165 


Overall  Ply  Orientation  Angle  (Alpha  +  Theta)  for  Top  and  Bottom 

Plies 


Figure  5.12:  Comparison  of  the  First  Mode  for  Varying  Thicknesses. 


Modal  Freq.  vs.  Overall  Ply  Orientation  Angle  2nd  Mode 


Overall  Ply  Orientation  Angle  (Alpha  +  Theta)  for  Top  and  Bottom  Plies 


Figure  5.13:  Comparison  of  the  Second  Mode  For  Varying  Thicknesses. 
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Table  5.7:  Effect  of  Laser  Cut  Angle  and  Thickness 


Effect  of  Laser  Cut  Angle  and  Thickness  on  Wing  Modal  Frequencies 


Laser  Cut  Angle 

a 

-6 

-5 

-4 

-3 

-2 

-1 

0 

1 

2 

3 

4 

5 

6 


1st  Modal  Freci. 

(Hz) 

135  micron  Thick 

46.1 

47.6 

49.2 

50.8 
52.5 

54.1 

55.7 
57.0 

57.9 
58.0 

57.4 

56.2 

54.8 


1st  Modal  Freci. 

(Hz) 

165  micron  Thick 

58.1 

60.0 

61.9 

63.9 

65.9 

67.9 
69.6 

70.8 

71.3 

70.9 

69.9 

64.4 

66.8 


2nd  Modal  Freq. 

(Hz) 

135  micron  Thick 

71.2 

70.8 

70.1 

69.2 

68.2 

67.2 
66.1 

65.2 

64.5 

64.3 

64.5 
65.0 

65.5 


2nd  Modal  Freq. 

(Hz) 

165  micron  Thick 

83.8 

83.5 
83.0 

82.4 

81.6 

80.8 
80.0 

79.5 

79.5 
80.0 

80.8 

81.5 

82.2 


Based  on  the  results,  it  can  be  seen  that  there  is  a  significant  difference  present 
when  comparing  the  both  the  first  and  second  modal  frequencies.  The  area  between  the 
two  curves  represents  possible  modal  frequencies  that  could  be  exhibited  in  the 
engineered  wing  based  on  the  current  composite,  and  its  variations  due  to  the 
manufacturing  process.  The  difference  seen  is  similar  to  difference  seen  in  Section  5.2.1, 
as  the  thickness  of  the  composite  increased,  the  modal  frequency  increased.  However,  it 
should  be  noted  that  as  the  thickness  of  the  composite  increases,  the  peak  seen  in  the 
curve  for  the  first  modal  frequency  and  the  low  point  of  the  second  modal  frequency 
curves  shifts  from  approximately  a  3°  a  for  the  135  micron  composite,  to  approximately 
2°  a  for  the  165  micron  sample. 

The  mean  modal  frequency  values  for  this  run  is  shown  in  Table  5.8,  and 
compared  to  the  ideal  (using  specific  composite  material  property  values)  values  and  the 
experimental  results.  The  results  show  that  given  the  composite  idiosyncrasies,  the  first 
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mean  first  mode  will  be  lower  than  the  ideal,  and  the  mean  second  mode  higher  than  the 
ideal.  This  follows  the  trend  determined  by  the  experimental  values. 

Table  5.8:  Comparison  of  Modal  Results. 


Modal  Frequencies 

First  Mode 

Second  Mode 

Difference 

First  Mode 

Second  Mode 

Result  Ideal  Experimental 

60.1  (Mean)  Hz  62.6  Hz  58.1  Hz 

74.3  (Mean)  Hz  73.0  Hz  80.3  Hz 

4.0%  -  7.2% 

1.8%  -  10.0% 

5.4.  Matching  the  FEA  to  Experimental  Results 

With  such  a  wide  range  of  possible  values  for  the  modal  frequencies  of  the  wing, 

it  is  highly  probable  that  the  FEA  model  could  be  made  to  match  the  experimental  results. 
Figure  5.14  shows  the  results  from  Section  5.3  (Figures  5.12  and  5.13)  plotted  on  a  single 
graph.  Also  plotted  are  two  straight  lines  showing  the  first  and  second  modal  frequencies 
achieved  experimentally. 
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Modal  Freq.  vs.  Overall  Ply  Orientation  Angle 
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Figure  5.14:  Modal  Frequency  Results  vs.  Overall  Ply  Orientation. 

Based  on  these  results,  a  trial  and  error  method  was  used  in  an  attempt  to  match 
the  experimental  results.  Using  an  a  of  4.5°,  resulting  in  a  ply  orientation  of  [-4.5/85.5/- 
4.5],  and  thickness  of  158  microns,  a  first  modal  frequency  of  58.0  Hz,  and  a  second 
modal  frequency  of  80.3  Hz  was  achieve.  This  represents  a  0.17%  difference  between 
the  first  modal,  and  a  0%  difference  between  the  second  modal  frequencies  of  the 
experimental  compared  to  the  analytical  FEA  results.  Both  of  these  values  also  fell 
within  the  likely  limits  measured  within  the  carbon  fiber  samples.  In  addition  the  mode 
shapes  matched  the  same  flap  and  feather  modes  seen  experimentally,  Figures  5.15  and 


.Mode  1, 135  Microns 
.Mode  1, 165  microns 
.  Mode  2, 135  microns 
.Mode  2, 165  microns 
Experimental  First  Mode 
Experimental  Second  Mode 
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5.16.  This  shows  that  the  stiffness  characteristics  of  the  model  match  the  stiffness 


characteristics  of  the  experimental  analysis. 


Figure  5.15:  First  Mode  of  Wing  Matched  to  Experimental  Results. 
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Figure  5.16:  Second  Mode  of  Wing  Matched  to  Experimental  Results. 


5.5.  Further  Discussion 

Based  on  the  results  seen  with  the  FEA  analysis,  it  can  be  determined  that  the 
errors  caused  in  the  manufacturing  process  play  a  significant  role  in  detennining  the 
modal  response  of  the  engineered  wing  as  would  be  expected.  While  for  the  case  of  the 
beam  specimens  in  Chapter  3,  this  caused  a  significant  reduction,  34.0%,  in  the  modal 
frequency  of  the  beam  compared  to  the  ideal  case,  the  curvature  of  the  veins  resulted  in  a 
lower  overall  percentage  reduction  in  the  modal  frequencies,  up  to  10%,  although  there  is 
still  a  wide  range  of  expected  frequencies. 

If  a  better  manufacturing  process  could  be  developed,  less  variation  in  the  fiber 
angle,  and  the  thickness  of  the  composite  could  significantly  reduce  the  variation, 
although  due  to  the  small  tolerances  present  in  the  manufacturing  of  the  composite,  some 
variation  would  be  inevitable  using  the  hand  lay-up  technique  required  here  at  AFIT. 
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Another  method  to  control  the  fiber  orientation  would  be  to  choose  a  ply  orientation 
further  away  from  the  0°  that  is  present  in  two  plies  of  the  composite.  This  would  reduce 
the  stiffness  of  the  model,  but  as  seen  when  varying  the  90°  ply,  less  variation  would 
occur.  However  this  would  not  improve  the  overall  stiffness  of  the  wing  without  adding 
more  plies,  and  thus  more  mass. 

The  FEA  model  does  match  the  engineered  wing,  however  the  engineered  wing 
does  not  match  the  biological  wing,  Table  5.9.  This  indicates  that  further  design  work 
still  needs  to  be  done  on  the  wing.  Given  the  near  infinite  number  of  possible 
arrangements  of  the  composite  material,  possible  changes  in  vein  width  to  the  wings,  and 
the  addition  of  camber  to  wing,  the  potential  for  an  engineered  wing  made  out  of  the 
composite  is  not  necessarily  unobtainable.  This  is  one  of  the  main  reasons  that  the 
MATLAB  code  was  created  in  order  to  generate  the  FEA  model.  Now  that  the  FEA 
model  has  been  matched  experimentally,  using  the  knowledge  gained  based  on  the  effects 
of  the  composite  material,  it  can  be  easily  modified  and  further  iterations  on  the  design 
can  be  tested  prior  manufacturing. 


Table  5.9:  Examination  of  Results. 


Examination  of  Results 

First 

Mode 

Second 

Mode 

Difference  from 
Biological  1st  Mode 

Difference  from 
Biological  2nd  Mode 

Biological 

86  Hz 

106  Hz 

- 

- 

Experimental 

58.1  Hz 

80.3  Hz 

32.44% 

24.25% 

FEA  Min 

46.08  Hz 

64.26  Hz 

46.42% 

39.38% 

FEA  Max 

71.3  Hz 

83.8  Hz 

17.1% 

20.94% 
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VI.  Summary  and  Conclusions 


The  use  of  unidirectional  carbon  fiber  composite  material  in  the  manufacturing  of 
FMAV  wings  provides  a  lightweight  material  with  a  high  specific  strength.  Previous 
research  done  at  Harvard  University  on  FMAV  wings  has  shown  that  the  material  when 
laid  in  a  [0/90/0]  orientation  was  capable  of  producing  a  viable  wing  for  use  in  a  MAV. 
Since  metals  and  plastics  have  proven  to  either  have  insufficient  strength  or  to  have  too 
much  mass,  the  [0/90/0]  orientation  of  unidirectional  carbon  composite  was  chosen  to 
form  the  vein  structure  of  a  FMAV  wing  based  on  the  forewing  of  the  Manduca  Sexta. 

The  focus  of  this  research  was  to  develop  a  finite  element  model  capable  of 
accurately  predicting  the  observed  modes  of  an  engineered  Manduca  Sexta  forewing. 
Due  to  considerable  variation  in  material  properties  of  the  [0/90/0]  carbon  fiber 
composite  used  in  the  manufacturing  of  the  vein  structure  within  the  wing,  a  considerable 
effort  was  put  forth  in  studying  the  effects  of  the  off  specification  composite.  The  FEA 
model  of  the  wing  was  used  to  study  the  effects  that  the  material  variation  within  the 
composite  would  have  on  the  modal  response  of  the  engineered  wing.  A  summary  of  the 
research  perfonned  for  this  project  is  stated  below. 

6.1.  Summary 

The  manufacturing  of  the  composite  material  was  perfonned  at  AFIT.  Sheets  of 
unidirectional  carbon  fiber  pre-preg  were  arranged  in  into  the  [0/90/0]  orientation  using  a 
hand  lay-up  technique.  This  lay-up  was  then  cured  in  a  heatpress  in  order  to  cure  the 
epoxy  and  fonn  the  composite  material.  A  laser  was  then  used  to  cut  the  composite  into 
the  desired  vein  structure  shape.  Using  beam  test  specimens,  the  modal  frequency  of  the 
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composite  was  evaluated  using  a  SLV.  However,  the  results  of  the  SLV  tests  proved 
inconsistent  with  FEA  results  of  the  beam  specimens. 

Since  the  composite  material  has  been  verified  using  ASTM  D  3039  and  ASTM 
D  3518  standards  using  20  ply  test  specimens,  errors  in  the  manufacturing  process  were 
investigated  in  order  to  determine  the  inconsistent  results.  Using  experimental  methods, 
the  thickness,  ply  orientation  and  laser  cut  angle  were  examined.  Small  variations  in  the 
thickness,  ply  orientation  and  laser  cut  angle  were  each  individually  determined  to  affect 
the  modal  response  of  the  beam  specimen. 

Experimental  modal  analysis  was  perfonned  on  test  specimens  using  a  SLV  in 
order  to  ensure  material  properties  of  the  composite  material.  An  average  value  of 
144.9Hz  was  found  for  the  first  modal  frequency  for  the  first  bend  mode  of  test 
composite  beam  specimens.  FEA  results  had  predicted  a  first  modal  response  of  219.5 
Hz.  An  investigation  into  the  likely  cause  of  such  a  difference  was  made.  The  thickness 
of  the  composite  was  measured  and  found  to  vary  between  135  and  165  microns.  Taking 
the  thickness  into  account,  the  numerical  frequency  of  the  beam  would  vary  from  197.6 
to  234.1  Hz,  a  significant  amount,  but  not  enough  to  account  for  the  144.9  Hz  average 
discovered. 

The  next  thing  looked  at  was  hand  lay-up  manufacturing  process  of  the  composite 
material.  Test  specimens  were  measured  using  an  optical  microscope  in  order  to 
detennine  the  true  orientation  of  the  fibers.  It  was  found  that  the  placement  of  the 
specimens  in  the  laser  cutter  could  be  off  approximately  2.5°.  This  would  affect  the 
overall  ply  orientation  of  the  specimen.  Also  investigated  was  the  ability  of  the  hand  lay¬ 
up  process  to  accurately  manufacture  the  desired  ply  angle.  Investigation  of  test  samples 
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using  an  optical  microscope  again  found  an  approximately  2.5°  deviation  from  the 
specified  ply  angle. 

Taking  the  thickness,  laser  cut  angle  and  ply  orientation  of  the  laminate  into 
account,  a  Monte  Carlo  solution  was  performed  in  order  to  detennine  the  probable  and 
possible  first  modal  frequencies  of  the  beam.  The  Monte  Carlo  solution  detennined  that 
a  mean  value  of  154.2  Hz  was  the  most  probable  for  the  modal  frequency  of  the  beam, 
with  the  majority  of  solutions  varying  between  100  and  200  Hz.  The  144.9  Hz  average 
was  well  within  the  bell  curve  of  the  solution,  showing  that  there  was  indeed  off 
specification  designs  in  the  composite  material,  and  that  it  would  affect  the  modal 
frequency  of  structures  made  of  the  material. 

Since  the  manufacturing  of  the  composite  material  can  create  variations  in  the 
specified  material  properties  of  the  composite,  the  effects  that  the  variation  of  thickness, 
laser  cut  angle,  and  ply  orientation  were  examined  using  FEA  for  the  engineered  wing. 
Using  MATLAB,  input  files  were  created  for  ABAQUS  in  order  to  generate  multiple 
FEA  model. 

An  FEA  model  of  the  engineered  wing  was  constructed  using  a  developed 
MATLAB  code.  The  vein  structure  was  modeled  using  B32  elements  and  the  membrane 
using  S8R  and  STRI65  elements.  MATLAB  was  chosen  in  order  to  more  easily  render 
the  material  properties  of  the  beam  element  representing  the  composite  vein  structure  of 
the  wing,  and  in  order  to  quickly  and  easily  incorporate  the  composite  variations  into 
multiple  FEA  simulations.  Each  beam  in  the  model  possessed  a  unique  material  property 
and  profile  based  on  its  location  on  the  curvature  of  the  vein.  The  angle  of  the  element 
was  determined  in  the  global  axis  of  the  model,  and  then  the  material  properties  of  the 


120 


composite  material  was  determined  based  upon  that  angle  using  the  transformation 
equation.  A  homogeneous  material  property  for  the  beam  element  was  created  by 
changing  the  element  profile  and  density  using  the  effective  moment  of  inertia,  and  then 
calculating  the  mass  so  that  it  matched  that  of  the  original  element  with  the  physical 
square  cross  section.  In  addition,  since  the  current  engineered  wing  design  did  not  fully 
match  the  modal  results  of  the  biological  wing,  the  modal  was  made  so  that  the 
dimensions  of  the  veins  could  be  easily  modified,  as  well  as  the  material  properties. 

A  model  was  created  based  upon  the  specifications  for  the  engineered  wing. 
While  the  first  four  mode  shapes  matched  the  modes  determined  by  the  experimental 
analysis,  the  first  and  second  modal  frequencies  were  off  by  7.19%  and  10.0% 
respectively  from  the  58.2  and  80.3  Hz  measured  experimentally.  The  effects  of  off 
specification  composite  were  then  examined. 

As  with  the  composite  beam  specimen,  the  thickness,  laser  cut  angle  and  ply 
orientation  of  the  composite  material  was  varied  independently  in  order  to  determine  the 
effects  each  would  play  on  the  modal  response  of  the  wing.  The  thickness  was  varied 
from  135  to  165  microns,  resulting  in  a  modal  frequency  varying  between  55.6  and  69.6 
Hz  for  the  first  mode,  and  66.1  and  80.0  Hz  for  the  second  mode.  The  variation  of  the 
mid-ply  orientation  from  85°  to  95°  showed  less  than  1%  variation  in  the  frequency, 
while  variation  of  the  top  and  bottom  plies  orientation  from  -5°  to  5°  resulted  in  modal 
frequency  change  between  57.7  and  64.3  Hz  along  a  concave  down  curve  peaking  at 
approximately  2°  for  the  first  mode,  and  between  71.9  and  77.2  Hz  along  a  concave  up 
curve  with  a  valley  at  approximately  2°.  Varying  the  laser  cut  angle  resulted  in  similar 
results  to  those  varying  the  top  and  bottom  plies  orientation.  This  led  to  the  conclusion 
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that  the  effect  of  the  top  and  bottom  plies  orientation  and  thickness  of  the  composite 
played  the  most  significant  roles  in  detennine  the  modal  response  of  the  wing. 

The  laser  cut  angle  was  varied  between  -6°  to  6°  for  thicknesses  of  135  and  165 
microns.  This  resulted  in  a  curve  representing  a  range  of  probable  values  for  modal 
frequencies  that  could  be  exhibited  by  the  wing  based  on  observed  variations  in  the 
composite  material,  figure  5.13.  It  was  shown  that  the  experimentally  observed  modal 
frequencies  of  the  wing  fell  within  the  probable  range. 

A  trial  and  error  mode  was  used  varying  the  thickness  and  laser  cut  angle  of  the 
wing  in  order  to  match  the  experimentally  observed  modal  characteristics  of  the 
engineered  wing.  It  was  found  that  a  thickness  of  158  microns  and  a  laser  cut  angle  of 
-4.5°  resulted  in  a  modal  frequency  of  58.0  Hz  and  80.3  Hz  for  the  first  and  second  mode. 
This  represented  a  0.17%  and  0.0%  difference  from  the  experimental  results.  The  mode 
shapes  for  the  modal  also  matched  those  determined  experimentally. 

It  was  detennined  that  the  experimental  results  of  the  wing  fell  within  the  bounds 
established  by  the  variation  in  the  composite  material.  Using  a  trial  and  error  method,  an 
FEA  model  of  the  wing  with  a  +8  micron  thickness  and  a  -4.5°  laser  cut  angle  was 
created  that  matched  the  experimental  modal  results  of  the  engineered  within  1% 
variation,  thus  validating  the  assumption  that  the  manufacturing  of  a  [0/90/0]  orientation 
led  to  major  variations  of  the  wing  response. 
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6.2.  Conclusions 


6.2.1.  Carbon  Fiber  Material 

The  use  of  the  carbon  fiber  composite  resulted  in  a  wide  range  of  variability  in  the 
engineered  wing  predicted  by  the  FEA.  It  does  however  posses  a  high  specific  modulus 
not  available  in  any  isotropic  material.  This  allows  the  material  to  be  lightweight,  but 
still  have  a  high  strength.  Due  to  the  low  mass  required  for  the  wing,  composites  seemed 
to  be  the  only  viable  solution  due  to  the  inertial  loading  caused  by  the  flapping  of  the 
wing  [2],  However  this  does  not  mean  that  the  YSH-70-A/RS-3C  composite  currently 
used  is  the  best  solution. 

There  are  numerous  amounts  of  other  composite  materials  commercially 
available,  each  with  their  own  unique  material  property.  Likewise,  the  ply  lay-up  and 
number  of  plies  have  a  great  affect  on  material  properties.  The  [0/90/0]  ply  layup 
currently  used  has  proven  to  be  sensitive  to  ply  orientation  and  thickness  variations. 
Repeatable  production  of  an  engineered  wing  with  consistent  dynamic  properties  would 
not  be  possible  given  the  current  material  and  manufacturing  technique. 

However,  if  a  quasi-isotropic  material  could  be  found  that  would  still  meet  the 
mass  and  strength  criteria  required  by  the  MAV  design,  it  could  alleviate  some  of  the 
issued  faced  with  the  ply  orientation.  Quasi-isotropic  materials  are  unidirectional 
composite  materials  where  the  laminate  arranged  such  that  the  material  properties  in  x,  y, 
and  z  direction  are  the  same.  If  a  material  and  orientation  could  be  found  or  developed 
that  would  provide  the  same  axial  stiffness  (ability  of  a  beam  to  resist  bending)  in  both 
the  x  and  y  directions  of  the  composite,  then  a  composite  material  could  be  viable.  This 
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material  would  be  less  sensitive  to  variation  in  ply  orientation,  however  would  most 
likely  face  the  same  issues  due  to  variations  in  thickness. 

6.2.2.  Carbon  Fiber  Manufacturing  Process 

The  current  manufacturing  process  of  the  composite  material  has  been  attributed 
to  the  source  of  the  ply  orientation  error.  It  is  apparent  that  if  the  current  3-ply  [0/90/0] 
orientation  is  to  be  used,  the  manufacturing  process  needs  to  be  improved.  Given  the 
significant  variations  in  the  modal  response  of  the  wing  due  to  small  changes  in  angle,  in 
order  for  a  reliable  and  repeatable  process  for  manufacturing  the  wing,  a  tolerance  of  only 
1°  would  be  acceptable  for  the  total  ply  orientation  alignment  (laser  cut  angle  +  ply  layup 
misalignment).  Any  further  variation  of  the  ply  orientation  of  the  material  would  create  a 
variation  in  the  modal  response  beyond  2.5%.  This  seems  highly  unlikely  given  the 
experimental  results  based  on  the  hand  layup,  and  the  fact  that  this  does  not  consider  the 
variation  due  to  the  thickness  of  the  composite.  It  should  be  noted  that  a  variation  of  less 
than  5  microns  would  be  needed  to  maintain  the  same  2.5%  variation  in  modal  response 
frequency,  even  with  an  ideal  ply  orientation.  Given  the  large  variation  experimentally 
measured,  the  probability  of  creating  a  composite  that  would  have  a  small  variation  in 
thickness  and  ply  orientation  using  the  current  method  would  be  extremely  improbable. 
Instead,  some  type  of  automated  manufacturing  process  would  be  required  in  order  to 
achieve  such  tolerances. 
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6.2.3.  Monte  Carlo  Solution 


A  Monte  Carlo  solution  was  used  in  order  to  evaluate  the  wire  beam.  One  of  the 
reasons  behind  the  Monte  Carlo  solution  is  not  only  to  find  out  what  is  possible,  but  also 
what  is  probable.  While  a  Monte  Carlo  solution  does  provide  useful  information  as  to  the 
probability  of  results  expected  in  the  beam  specimen,  the  large  amount  of  data  required  in 
order  to  gain  a  significant  statistical  average  made  it  unfeasible  for  use  on  the  engineered 
wing. 

6.2.4.  Effects  of  Composite  Variation  on  the  Wing 

It  is  without  doubt  that  the  variation  in  the  composite  has  a  significant  effect  on 
the  modal  response  of  the  wing.  Based  on  such  a  wide  range  of  variation  that  was 
determined,  it  would  be  difficult  for  a  repeatable  wing  to  be  developed  using  the  current 
method  of  manufacturing.  However,  it  should  be  noted  that  the  mode  shapes  of  the  wing 
are  consistent,  meaning  that  the  response  of  the  wing  is  the  same,  just  the  numerical 
frequency  at  which  the  response  occurs  would  change.  This  means  that  while  the  wing 
does  not  match  the  stiffness  of  the  biological  wing,  it  shares  the  same  characteristics, 
meaning  that  the  wing  should  behave  similarly  to  the  biological  wing,  albeit  at  different 
frequencies.  Getting  the  wing  to  match  the  biological  wing  would  be  especially  be 
difficult  given  the  fact  that  while  the  wing  is  a  single  entity,  it  would  always  be  required 
to  work  in  pairs. 

This  would  add  a  complexity  of  having  to  “pair  up”  sets  of  wings  together  that 
share  similar  modal  response.  Failure  to  do  so  could  cause  significant  control  issues  for  a 
FMAV,  such  as  having  to  flap  two  wings  at  different  frequencies  in  order  to  generate 
equal  amounts  of  lift.  One  possible  solution  to  this  issue  would  be  to  use  vein  structures 
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that  are  cut  from  the  same  general  area  of  the  composite  material.  Since  the  fibers  are 
imbedded  in  the  pre-preg  of  the  epoxy  matrix  mechanically,  the  fibers  are  generally 
parallel  to  one  another.  Wings  that  are  cut  out  of  the  same  piece  of  composite  should 
have  less  variation  in  the  material  properties  compared  to  one  another  versus  those  cut 
from  other  sheets  of  the  composite.  This  can  be  evident  by  looking  at  the  beam  samples, 
which  were  cut  out  of  a  single  sheet  of  composite. 

However,  this  would  not  account  for  the  variation  in  properties,  only  amongst  the 
variation  between  wings.  It  can  be  concluded  that  this  would  only  be  a  feasible 
alternative  if  the  variation  in  the  MR  was  not  as  significant  as  long  as  matching  sets  of 
wings  shared  the  same  MR. 

6.2.5.  MATLAB  to  Generate  FEA  Model 

The  use  of  beam  elements  were  chosen  because  of  the  ability  for  the  geometry  to 
easily  be  modified  .  This  method  was  chosen  as  opposed  to  a  composite  shell  method 
where  the  number  and  orientation  of  the  plies  can  easily  be  changed,  but  the  geometry 
would  be  difficult  to  do  so.  This  decision  was  made  at  the  beginning  of  the  project  prior 
to  the  composite  material  variations  being  known,  when  it  was  thought  that  the  geometry 
of  the  wing  would  need  to  be  easily  changed  in  order  to  aide  in  the  design  process.  This 
would  allow  the  user  to  vary  the  width  of  the  veins  in  a  similar  matter  that  the  material 
orientation  was  varied.  Multiple  runs  with  various  designs  vein  widths  could  easily  be 
submitted  and  solved  to  determine  the  effect  that  a  change  in  vein  width  would  have  on 
the  overall  model.  While  this  is  still  a  useful  feature,  its  value  to  this  research  was 
limited  to  only  fine  tuning  the  model. 
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The  beam  elements  proved  to  be  faster  to  solve  than  the  shell  elements.  This  was 
especially  useful  for  cutting  down  the  time  it  took  to  solve  the  Monte  Carlo  solution  for 
the  sample  beam  specimens.  They  were  also  easier  to  generate  the  nodes  because  only  a 
single  line  was  required  as  a  reference  for  the  location  of  the  vein,  as  opposed  to  the  two 
that  would  be  required  for  shell  elements.  This  made  it  significantly  easier  if  the 
geometry  of  the  wing  was  required  to  change.  However,  the  downside  to  beam  elements 
was  that  each  beam  element  required  its  own  profile  and  material  property.  This  meant 
that  a  MATLAB  code  was  required  in  order  to  develop  the  node  sets,  element  sets, 
material  properties  and  bream  profiles.  This  meant  that  more  time  needed  to  be  spend 
“upfront”  in  developing  the  MATLAB  code,  as  opposed  to  using  the  CAD  program  in 
ABAQUS. 

The  MATLAB  code  made  the  wing  significantly  easier  to  modify,  enabling  the 
geometry,  beam  profile  and  composite  material  variations  to  easily  be  changed.  This 
means  that  the  code  can  easily  be  used  to  modify  the  width  of  any  of  the  veins,  or  even  be 
easily  modified  to  have  a  camber.  The  downside  to  the  present  code  is  that  due  to  the  I- 
beam  shape  of  the  composite,  the  beam  element  can  only  handle  3-ply  lamina.  This 
places  a  significant  limitation  on  the  code  and  its  capabilities,  considering  the  value  of  2 
or  4  plies  may  have  in  producing  a  more  quasi-isotropic  material.  This  feature  can  be 
changed  to  support  other  numbers  of  plies. 

6.3.  Suggested  Future  Work 

Based  on  the  results  of  the  analysis,  further  modification  of  the  engineered  wing  is 
required  in  order  to  match  the  characteristics  of  the  Manduca  Sexta  Forewing.  One  thing 
not  yet  examined  is  different  ply  angles.  The  equations  in  Chapter  2  show  that  the 
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material  properties  of  the  composite  change  significantly  as  the  ply  orientation 
approaches  0°.  Also,  the  results  from  Section  5.3,  also  confirm  that  the  modal  response 
of  the  wing  can  be  significantly  altered  based  on  the  ply  orientation.  Countless  other  ply 
orientations  could  possibly  be  used  in  the  manufacturing  of  the  vein  structure  that  would 
be  less  susceptible  to  the  idiosyncrasies  examined  for  this  project. 

Another  thing  that  can  be  examined  is  camber  of  the  wing.  Currently  the  wing  is 
created  on  a  flat  surface.  As  Sims  has  shown,  the  effect  of  camber  on  the  wing  will  have 
a  stiffening  effect  without  adding  any  additional  mass  [15].  Since  the  model  was 
developed  using  a  MATLAB  code,  imparting  a  camber  to  the  wing  could  feasibly  be 
done  relatively  easily. 

Focusing  on  possible  numerical  research,  a  design  should  be  created  that 
incorporates  the  camber  observed  on  the  biological  wing.  A  parametric  study  on  various 
ply  orientations  of  the  composite  material  and  vein  widths  should  be  conducted.  This 
could  be  easily  done  by  modifying  the  current  MATLAB  code  that  generated  the  FEA 
model  of  the  wing.  At  this  point,  it  is  clear  that  the  current  manufacturing  process  of  the 
composite  created  a  significant  variation  in  the  modal  response  of  the  wing.  A  less 
susceptible  design  needs  to  be  created,  whether  by  introducing  a  different  ply  orientation, 
a  camber,  or  using  a  different  composite  manufacturing  process  -  or  any  combination  of 
the  three,  in  order  to  create  a  viable  engineered  wing  that  can  be  optimized  to  matching 
the  biological  Manduca  Sexta  forewing  characteristics. 
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Appendix  A.  YSH-70-A/RS-3C  Material  Properties 
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NGF  High  Stiffness  Fibers 
Nominal  Unid  Laminate  Properties 


Appendix  B.  YSH-70-A/RS-3C  Material  Property  Calculation  Code 

This  code  is  the  MATLAB  code  that  was  used  to  solve  the  material  properties  of 
the  YSH-70-A/RS-3C  composite. 


function 

[wr, Ex, Ey, Exb, Eyb, Vxy, Vyx, Gxy, Gxz , Gyz , Gxyb, El , E2 , G12 ] =LAMINATE_NASA_HAL 
PIN_HYBRID_BASELINE_FUNC (l,b, h, alpha, theta_d,  prnt) 

o, 

o 

[wr, Ex, Ey, Exb, Eyb, Vxy, Vyx, Gxy, Gxz , Gyz , Gxyb, El ,  E2  ,  G12 ] =LAMINATE_NASA_HAL 
PIN_HYBRID_BASELINE_FUNC (40E-3, 2 . 5E-3,  140E-6,  0,  [0  90  0 ]  ,  1 ) 

%  Inputs: 

%  (1)  beam  length  (m) 

%  (b)  width  of  beam  (m) 

%  (h)  laminate  thickness  in  meters  (m) 

%  (alpha)  angle  of  laser  cut 

%  (theta_d)  laminate  layup  in  degrees  [0  90  0] , [0  45  -45  0) 

%  (prnt)  (0  =  no)  (1  =  yes)  to  print  values 


Outputs : 

[wr]  First  Resonant  Frequency  of  the  beam 

[Ex]  Laminate  Axial  Stiffness 

[Ey]  Laminate  90  Axial  Stiffness 

[Vxy]  Laminate  Poissons  Ratio 

[Vyx]  Laminate  Poissons  Ratio 

[Gxy]  Laminate  Shear  Modulus  x-y  direction 

[Gxz]  Laminate  Shear  Modulus  x-z  direction 

[Gxy]  Laminate  Shear  Modulus  y-z  direction 

[El]  Lamina  Axial  Modulus 

[E2]  Lamina  90  Axial  Modulus 

[G12]  Lamina  Shear  Modulus 


Input  Fiber  and  Matrix  Values 

(pa) 

(pa) 

(pa) 

(pa) 


Ef  =  6.55E11; 

O, 

o 

%Ef  =  9.23E11; 

g, 

o 

Em  =  2E9; 

Q, 

O 

Gf  =  6.57E9; 

g. 

o 

(Estimate) 

Gm  =  3.316E9; 

o, 

o 

3C)  -  http : //www31 . 

.  ocn .  n< 

vf  =  0.3; 

o, 

o 

vm  =  0.3; 

g, 

o 

(RS-3C) 

df  =  2140; 

g, 

o 

dm  =  1193; 

g, 

o 

Vf  =  0.639; 

g, 

o 

(estimate) 

Vm  =  1  -  Vf; 

g, 

o 

tf  =  7  E -  6 ; 

g, 

o 

Elastic  Modulus  of  Fiber  (YSH-70A) 
Elastic  Modulus  of  Fiber  (YSH-95A) 
Elastic  Modulus  of  Matrix  (RS-3C) 
*Shear  Modulus  of  the  Fiber 


(pa)  -  Shear  Modulus  of  the  Matrix  (RS- 

. jp/~ngf / english/product /pi . htm 

-  Poisson's  Ratio  of  the  Fiber  (YSH-70A) 

-  Poisson's  Ratio  of  the  Matrix 

(kg/mA3)  -  Density  of  the  Fiber  (YSH-70A) 
(kg/mA3)  -  Density  of  the  Matrix  (RS-3C) 

-  *Fiber  Volume  Fraction  (YSH-70A) 


(m) 

theta_d  +  alpha; 


d  = 

r  =  theta_d  *  pi  /  180; 


theta 
theta_ 
t=h/length (theta_r) ; 
Input 

H  =  h/2 ; 


-  Matrix  Volume  Fraction 

-  Fiber  Diameters 

Fiber  Orientation  +/-  Laser  Cut  angle 

-  Layer  Fiber  Orientation  in  Radians 
Lamina  Thickness  based  on  Laminate  Thickness 


Half  laminate  thickness  (Fig  5.2  Herakovich) 
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NL  =  length (theta  r) ;  %  Number  of  Layers 

Le  =  b./sin(theta  r) ;  %  (m)  -  Effective  Fiber  Length  based  on  theta 
and  beam  width 

%%  Make  Sure  Le  is  calculated  properly 
for  x=l : NL 

if  theta  d(x)  ==  0,  LE  =  1;  end  %  -  Account  for  Divide  by  0 

if  Le(x)  >  l,Le(x)=l;  end  %  -  This  is  ~0-4  degrees  @  2.5  mm 

end 

%%  Halpin-Tsai  Calculations  -  Account  for  Fiber  Orientation  and 
Effective  Fiber  Length 

%  pg  350  first  paragraph 

http : // www . abdmatrix . com/ phcdl/ upload/ f undamentals/Halpin- 
Tsai%20Equations-A%20Review . pdf 

%  The  hybrid  implmentation  of  this  only  uses  El  from  Halpin-Tsai, 
other 

%  values  follow  the  standard  calcs  from  Herakovich 
%  the  values  suggested  by  the  Halpin-Tsai  paper  produced  E2,G12 
values 

%  that  were  excessively  high  based  on  the  estimate  of  [2  1]  for 
zeta . 

MAT=zeros (NL, 6) ;  %  Pre-Allocate  MAT  ARRAY 

if  prnt==l 

fprintf ( 'Halpin-Tsai  Calculations'^  1 ) 

end 

for  x=l : NL 

zeta (x, : )  =  [2*Le(x)/tf  1E6  1E6] ;  %  Halpin-Tsai  zeta  factors 
[zetal  zeta2  zeta3] 

ada (x,  : )  =  (Ef-Em)  . / (Ef tzeta (x,  : ) *Em)  ;  %  ada=[adal  ada2  ada3] 

E=Em. * (1+zeta (x) . *ada (x) *Vf ) . / (1-ada (x, : ) *Vf ) ;  %  E=[E1  E2  G12] 
E1=E(1);  %  Only  use  El,  ignore  E2  and  G12 

E2=Em/ (Vf * (Em/Ef-1 ) +1 ) ;  %  ( 1 1 . 1 4  -  Herakovich) 
vl2=Vf* (vf-vm) +vm;  %(11.9  -  Herakovich) 

v21=E2/El*vl2;  %pg4  -  http://pas.ce.wsu.edu/CE537- 

1 / Lectures /Rule%2 Oof %20Mixtures .pdf 

G12  =  l/ (Vf/Gf+Vm/Gm) ;  %  (11.23  -  Herakovich) 

G13=G12 ; 

G23=E2/2* ( l+vl2 ) ; 

Den=df *Vf tdm*Vm; 

MAT (x, : ) = [ E 1  E2  G12  vl2  v21  Den];%Store  Lamina 
if  prr.t—  1 

fprintf  ('  E1=%E  E2=%E  vl2=%f  G12=%E  Den  =  %E\n’,... 

[El  E2  vl2  G12  Den] ) 

end 

end 

%%  Compute  Layer  Height  Data 
range=zeros (1,NL+1) ; 
for  x  =  1 : NL+1 
if  x  ==  1 

range (x)  =  NL*t/2; 

else 

range (x)  =  NL*t/2-t* (x-1)  ; 

end 

end 
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%%  Compute  Qp  Matrix 

%  Q^Global  based  on  varying  Zeta  Values  due  to  Le  differences  for 
each  Lamina 

Q= zeros ( 3 , 3 , NL) ; 
for  x  =  1 : NL 

E 1 =MAT (x, 1) ; E2=MAT (x, 2) ;G12=MAT (x, 3) ;vl2=MAT (x, 4) ;v21=MAT (x, 5) 
Qll  =  El/ (1-v12*v21) ;  %(4.16  -  Herakovich) 

Q12  =  (vl2*E2) / (1-v12*v21) ;  %(4.16  -  Herakovich) 

Q22  =  E2/ (1-v12*v21) ;  %(4.16  -  Herakovich) 

Q66  =  G12;  %(4.16  -  Herakovich) 

Q ( :  ,  :  ,  x)  =  [Qll  Q12  0;Q12  Q22  0;0  0  Q66];%Global  Q  for  Lamina 

end 

%  Init  A, B,D  matrices 
A=zeros (3,3) ; 

B=zeros (3,3) ; 

D=zeros (3,3) ; 

%  Q_ply  and  A, B,D  Matrices 
for  x=l : NL 

m=cos (theta_r (x) ) ; 
n=sin (theta_r (x) ) ; 

T= [mA2  nA2  2*m*n;nA2  mA2  -2*m*n;-m*n  m*n  mA2-nA2];%  Motavalli 

pl37 

QP=  T*Q ( : , : , x) *T ' ;  %  Motavalli  pl37  -  QP  is  different  for  each 

lamina 

A=A+QP* (range (x) -range (x+1) ) ; 

B=B+1/2*QP* (range (x) A2-range (x+1) A2) ; 

D=D+1/3*QP* (range (x) A3-range (x+1) A3) ; 

end 


%%  Compute  Global  Laminate  Properties 


%  h  |  | - | - >  y  (b  =  width) 

%  | _ |  (h  =  height) 

%  b 

%  X-axis  is  along  the  length  of  the  beam 
%  Y-Axis  is  along  the  width  of  the  beam 
%  Z-Axis  is  along  the  thickness  of  the  beam 

%%  Compute  Global  Laminate  Properties  -  Axial  (NASA) 

ABD= [A  B;B  D] ;  %  Assembled  ABD  Matrix 

Ex=l/h*det (ABD) /det (ABD ( [2  3  4  5  6], [2  3  4  5  6]));  %  Simplified 
Using  Cofactor  Expansion  about  1  -  NASA(84) 

Ey=l /h*det (ABD) /det (ABD ( [ 1  3:6], [1  3:6]));  %  Simplified  Using 
Cofactor  Expansion  about  2  -  NASA(85) 

astar=2*H*inv (A) ;  % ( 5 . 8 0  -  Herakovich) 

Vyx=-astar (1, 2) /astar (1, 1) ;  %  Motavalli  plOO 
Vxy=-astar (1 , 2) /astar (1 , 1) ;  %  Motavalli  plOO 

Gxy=l /h*det (ABD) /det (ABD ( [ 1  2  4  5  6 ] , [ 1  2  4  5  6]));  %  Simplified 
Using  Cofactor  Expansion  about  3  -  NASA(89) 

Gxz=Gxy; 

Gyz=Ey/2* (1+Vxy) ; 
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%%  Compute  Global  Laminate  Properties  -  Bending  (NASA) 
Exb=12/hA3*det (ABD) /det (ABD ( [1  2  3  5  6],[1  2  3  5  6])); 
Simplified  Using  Cofactor  Expansion  about  4 

Eyb=12/hA3*det (ABD) /det (ABD ( [1  2  3  4  6],[1  2  3  4  6])); 
Simplified  Using  Cofactor  Expansion  about  5 

Gxyb=12/hA3*det (ABD) /det (ABD ( [1  2  3  4  5]  ,  [1  2  3  4  5])) 
Simplified  Using  Cofactor  Expansion  about  6 
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Appendix  C.  Beam  .inp  File  Generator  Function 

This  Code  was  developed  to  generate  the  .inp  fde  for  the  wire  beam  used  in  Chapter  3. 

function  [Eigenvalue]  =Vary  Properties^Beam (alpha,  theta  original,  h,  t) 
%Frequency 

num  modes  =  1;  %the  number  of  modes  you  wish  to  find 
clamped  nodes  =  1;  %nodes  you  wish  to  clamp  i.e.  the  base 
hl=  h; 

den_cf  =  1750;%  density  of  the  carbon-fiber  composite 
%%  Variables 
%scale  nodes 

%note,  made  root  of  wing  =  0,0. 
x_offset=0;  %move  nodes  in  x  dir 
y^of f set=0 . 0 ;  %move  nodes  in  y  dir 

wing  length  =  0.04;  %desired  length  of  wing  in  meters 
%file  name 

File  =  'Beam  Test  out';  %  Name  of  the  file  created,  remember  to  change 
%  on  final  line  of  code  for  evaluation 

Node  File  =  'Beam_Test';  %  Name  of  file  where  nodes  are  located 
NodesMax  =  0;  %number  of  nodes  that  exist  (claculated  later) 

%%  Open  the  File 

fid2=fopen ( [File  ' . inp ' ] , ' w ' ) ;  %  Open  file  for  writing,  new  file 
%%  Heading 

%  Standard  Abaqus  Heading 
fprintf ( f id2  ,  ' * Heading \n ' ) ; 

fprintf ( f id2 , '**  Job  name:  BEAM  MODAL  BASELINE  Model  name:  %s\n',date) 
%%  Contact 

%  Contact  within  model 

fprintf (fid2,  ' *Preprint,  echo=NO,  model=NO,  history=NO,  contact=NO\n ' ) 
fprintf ( f id2 , ' **\n ' ) ; 

%%  Parts 

%  Parts  (doesn't  change  for  this  case) 
fprintf (fid2,  ' **  PART S \ n ' ) ; 

fprintf (fid2, ' *Part,  name=WING  FIT  SHOULDER  EE  %s\n',date);  %  name 

called  up  in  instance 

fprintf ( f id2 , ' **\n ' ) ; 

fprintf ( fid2 ,' *End  Partin'); 

fprintf ( f id2 , ' **\n ' ) ; 

fprintf (fid2, '**'); 

%%  Assembly 

%insert  the  Assembly  section  (doesn't  change  for  this  case) 
fprintf (fid2, ' **  ASSEMBLY\n ' ) ; 
fprintf ( f id2 , '**\n'); 

fprintf ( fid2 , ' *Assembly,  name=Assembly\n ' ) ; 
fprintf ( f id2 , ' **\n ' ) ; 

%%  Instance 

%  Insert  the  Instance  section 
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fprintf (fid2, ' *Instance, name=Part-l-l , 

part=WING  FIT  SHOULDER  EE  %s\n '  ,  date) ; %  name  called  from  parts 

%%  Insert  the  nodes  from  the  veins 
%Standard  Heading 

fprintf ( fid2 *Node ') ;  %heading  for  the  nodes 
fprintf ( f id2 , ' \n ' ) ; 

%  Take  the  nodes  out  of  another  inp  file. 

f idl=f open ( [Node^File  '  . inp ' ] ,  ' r+  '  )  ;  %  Open  file  for  reading/writing 
idx  =  0 ; 

while  ~feof(fidl) 
idx  =  idx  +  1; 
tline=fgetl (fidl) ; 

if  isempty (strfind (tline, ' *Node ' ) )  ==  0  %  Scan  till  *Nodes 

for  i=l : 2000000000 

tline=fgetl (fidl) ; 

if  isempty (strfind (tline, '*') )  ==  0%  Scan  till  *  Found 

NodesMax=i-l ; %calculate  the  new  value  of  maximum  nodes 
break 

else 

nodes (i, : ) =sscanf (tline,  %d,%f,%f,%f') 

end 

end 

break 

end 

end 

fclose (fidl ) ; 

%  for  loop  to  print  nodes  into  the  .inp  file 
%  nodes  will  be  scaled  at  this  time 

node^scale  =  wing_length/ (max (nodes ( : , 2 ) ) ) ;  %sets  scale  =  to  desired 
legnth 

fprintf ( fid2 ,' %d,  %2.10f,  %2.10f,  %2 . 10f\n  , [nodes (:, 1 ),.. . 

(nodes ( :  ,  2 ) -x^of f set) *node_scale, 

(nodes ( : , 3) -y_offset) *node_scale, 
nodes ( :  ,  4 )  ]  '  )  ; 

%%  Elements 
%  For  beam  elements 

fprintf (fid2, ' *Element,  type=B32 ' ) ;  %heading  for  the  elements 
fprintf ( f id2 , ' \n ' ) ; 

%Search  Input  File  for  *Elem  Locations 

fidl=fopen ( [Node  File  ' . inp ' ] , 1 r+ ' ) ;  %  Open  file  for  reading/writing 
idx  =  0 ; 

while  ~feof(fidl) 
idx  =  idx  +  1; 
tline=fgetl (fidl) ; 

if  isempty (strfind (tline, ' *Element,  type=B32 ' ) )  ==  0  %  ' ' 
for  i=l:NodesMax  %  Scan  till  *  Found 

tline=fgetl (fidl) ; 
if  isempty (strfind (tline,  '*') )  ==  0 

break 

else 
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elems . B32 (i,  :)=sscanf(tline,  '%f,%f,%f,%f,%f,%f,%f,%f,%f') 
end 

end 

break 

end 

end 

fclose ( f idl ) ; 

%  Write  nodes  into  the  program  file 

fprintf (fid2, ' %d, %6 . Of , %6 . Of , %6 . Of \n '  , [elems. B32] ' ) ; 

%%  Read  Stringer  Node  Sets 

%Search  Input  File  for  Stringer  Nsets 

fidl=f open ( [Node  File  ' . inp ' ] , ' r+ ' ) ;  %  Open  file  for  reading/writing 
for  j=l:12  % (Number  of  Stringers) 
idx  =  0; 

Stringers ( j ) . nset= [ ] ; 
while  ~feof(fidl) 
idx  =  idx  +  1; 
tline=fgetl (fidl) ; 

if  isempty (strfind (tline, ['*Nset,  nset=  Stringer-'  num2str(j) 
'  ,  internal ']))==  0  %  '  ' 

for  i=l:NodesMax  %  Scan  till  *  Found 

tline=fgetl (fidl) ; 

if  isempty (strfind (tline,  '*') )  ==  0 
break 

else 

Stringers (j) . nset= [ Stringers ( j ) .nset 
sscanf (tline, ' %f , ' ) ' ] ; 

end 

end 

break 

end 

end 

end 

fclose (fidl ) ; 

%%  Calculate  Angle  between  the  Stringer  B32  Elements 
%  needed  to  determine  material  properties, 
for  m  =  1 : ( length (elems . B32 ( : , 1 ) ) ) 

%determine  the  first  node  that  is  part  of  the  element  i.e  x 
coordinate 

node^str  =  elems . B32 (m, 2 ) ; 

%determine  the  second  (mid  point)  node  of  the  element  i.e  y 
coodiinate 

node  end  =  elems . B32 (m, 3 ) ; 

%location  of  first  node 
loc_str  = 

[ (nodes (node_str , 2 ) ) , (nodes (node_str , 3 ) ) , (nodes (node_str , 4 ) ) ] ; 

%  location  of  the  end  node 
loc_end  = 

[ (nodes (node_end, 2 ) ) , (nodes (node_end, 3 ) ) , (nodes (node_end, 4 ) ) ] ; 

%  determine  the  slope  of  the  line 

slope_2  =  ( loc  str ( 1 , 2 )  -  loc  end ( 1 , 2 ) ) / ( loc_str ( 1 , 1 )  - 

loc  end (1,1) ) ; 
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%determine  the  angle  with  horizontal  [inf, 0,0]  vector  from  origin 
theta_elem  =  atand ( slope_2 ) ; 

%store  it! 

angle (m, 1 ) =m;  %stores  the  row  number  in  angle  in  1st  column 
angle (m, 2) =  elems . B32 (m, 1 ) ;  %the  2nd  column  hold  element  number 
angle (m, 3 ) =theta  elem;  %  3rd  column  holds  the  angle 

end 


%%  Write  Nset  and  Elsets 
%  For  the  beam 

set_ctr=0;  %counts  the  number  of  Node  and  Element  Sets 
%  Create  a  new  section  for  each  Vein  element 
for  x=l : length (angle) 

%create  heading  and  pick  a  Nset 
fprintf (fid2, ' *Nset,  nset=  PickedSet%-2 . Of , 
internal \n ' , angle (x, 1 ) )  ; 

nset  nodes  =  sort (elems . B32 (x, 2 : 4 ) ) ; %sort  out  the  nodes  in 
asscending  order 

%print  nodes  in  that  Nset  i.e.  for  each  element 
fprintf (fid2, ; %6. Of,  %6.0f, 

%6 . Of \n ' , nset_nodes (1) , nset_nodes (2) , nset_nodes (3) ) ; 

%  for  Elsets,  create  elset  same  number  as  corresponding  Nset 
fprintf (fid2, ' *Elset,  elset=  PickedSet%-2 . Of , 
internal \n  , angle (x, 1 ) ) ; 

fprintf (fid2,  ; %6 . Of \n ' , angle (x, 2 ) ) ; 

set_ctr=set_ctr+l ;  %add  to  Set  counter 

end 


%%  Beam  Sections 

%  The  sections  will  start  with  the  veins  and  begin  with  Section  1  (B32) 

for  x=l : length (angle)  %  new  section  for  each  element 

fprintf ( fid2 ,'* *  Section:  Section-%-2 . Of  Profile:  Profile-%- 

2 . Of \n ' . . . 

,x,x);  %  create  section/profile 
fprintf (fid2, . . . 

' *Beam  Section,  elset=  PickedSet%-2 . Of ,  material=Material-%- 
2. Of,  temperature=GRADIENTS ,  section=I\n ' . . . 

,x,x);  %  match  the  element  with  it's  own  elset  and  material 

b (x) =2 . 5e-3;  %artif icial  width  of  the  beam  element 

%  Generate  the  I  beam  sections 

theta  =  angle (x, 3 ) +theta  original ( 1 ,:) ;  %determine  theta 
%For  the  top  ply  run  a  function  to  determine  Matl  Prop 
[ ~ , Ex, Ey, ~ , ~ , Vxy, Vyx, Gxy, Gxz , Gyz , ~ , ~ , ~ ,  ~ ]  =  .  .  . 

LAMINATE_NASA_HALPIN_HYBRID_BASELINE_FUNC (4  0E- 
3,  b (x) , h, alpha, theta (1,1) , 0 ) ; 

%  Capture  the  Mat  Properties  so  we  don't  have  to  run  this  again 
Matl  Ex(x,l)  =  Ex;  Matl  Ey(x,l)  =  Ey;  Matl  Vxy(x,l)  =  Vxy; 

Matl_Vyx (x, 1 )  =  Vyx; 

Matl  Gxy(x,l)  =  Gxy;  Matl_Gxz (x, 1 )  =  Gxz ; Matl^Gyz (x, 1 )  =  Gyz; 

%Second  ply 

[~ , Ex, Ey, ~ , ~ , Vxy, Vyx, Gxy, Gxz , Gyz , ~ , ~ , ~ ,  ~ ]  =  .  .  . 
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LAMI NATE_NAS A_HAL  P I N_H YBRI D_BASE L I NE_FUNC ( 4  0  E - 
3 , b (x) , h, alpha, theta (1,2) , 0 )  ; 

Matl  Ex(x,2)  =  Ex;  Matl  Ey(x,2)  =  Ey;  Matl  Vxy(x,2)  =  Vxy; 

Matl  Vyx(x,2)  =  Vyx; 

Matl_Gxy (x, 2 )  =  Gxy;  Matl_Gxz (x, 2 )  =  Gxz ; Matl_Gyz (x, 2 )  =  Gyz; 

%Third  ply 

[ ~ , Ex, Ey, ~ , ~ , Vxy, Vyx, Gxy, Gxz , Gyz , ~ , ~ , ~ , ~ ] = . . . 

LAMINATE_NASA_HALPIN_HYBRID_BASELINE_FUNC (4  0E- 
3, b (x) , h, alpha, theta (1,3)  ,  0 )  ; 

Matl  Ex(x,3)  =  Ex;  Matl  Ey(x,3)  =  Ey;  Matl  Vxy(x,3)  =  Vxy; 

Matl  Vyx(x,3)  =  Vyx; 

Matl_Gxy (x, 3 )  =  Gxy;  Matl_Gxz (x, 3 )  =  Gxz ; Matl_Gyz (x, 3 )  =  Gyz; 

%  Use  largest  value  of  Ex  to  decide  value  for  the  effective  modulus 
%  This  prevents  the  width  of  some  sections  from  being  too  large. 

Eff  Mod  =  max (Matl  Ex (x,  :));  %determines  max  Ex 

width (x,l)=  Matl_Ex (x, 1 ) /Ef f  Mod;  %makes  the  width  of  each  section 
a  scaler,  with  largest  being  =  1 

width  (x,  2)  =  Matl__Ex  (x,  2  ) /Ef  f  Mod; 
width  (x,  3)  =  Matl__Ex  (x,  3) /Ef  f  Mod; 

%  Prevent  I  (12) *2  from  being  >  I (11) +1(22) 
if  width (x,l)<  width (x, 2) /7 . 5 
width (x, 1 ) =width (x, 2 ) ; 

h  new=h/ (3- (width (x, 2 ) / (min (width (x,  :))))); 
h=h  new; 

end 

if  width(x,3)<  width (x, 2) /7 . 5 
width (x, 3 ) =width (x, 2 ) ; 

h  new=hl/ (3- (width (x, 2) / (min (width (x,  :))))); 
h=h  new; 

end 

%  Calculate  the  effective  density  for  the  Beam  Sections 
real  area  =  hl*b(x);  %real  cross  sectional  area  of  the  element 
faux_area  =  h*b (x) * (width (x, 1) /3  +  width (x, 2) /3  +  width (x, 3) /3) ; 
%area  of  I-beam  section 

area_scale  =  real_area/faux_area;  Iscale  factor  to  multiply  density 
by  so  mass  is  same  for  all  elements 
Density (x) =den_cf *area_scale; 

%  prints  out  in  (1,  h,  bl (base  width) ,  b2  (top  width) ,  tl  (base) , 
t2 (top) ,  t3 (mid  width) ) 

fprintf ( f id2 ,  '%  #6.12g,  %  #6.12g,  %  #6.12g,  %  #6.12g,  %  #6.12g,  % 

#6 . 12g,  %  #6 . 12g\n ' , . . . 

h/2 , h, b (x) * width (x, 3 ) , b (x) * width (x,l) ,h/3,h/3,b(x) * width (x,  2 ) )  ; 
%prints  out  the  orientation  of  the  beam,  same  for  all  cases 
fprintf ( fid2 ,' 0 ., 1 ., 0 .  \n'); 

%  put  thickness  back  to  original  thickness 
h=hl  ; 

end 


%%  End  Instance 

fprintf (fid2, ' *End  InstanceVn' ) ; 
fprintf ( f id2 , ' **\n ' ) ; 
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%%  Nset  for  boundry  conditions 

%create  a  set  of  nodes  where  the  boundry  conditions  will  be  applied 
fprintf (f id2,  ' *Nset,  nset=  PickedSet%-2 . Of ,  internal,  instance=Part- 1- 

1 \n '  .  .  . 

,  (set_ctr+3) )  ; 

%Put  in  the  nodes  you  wish  to  clamp  in  this  case  the  base 
fprintf (fid2, ' %d, \n ' , clamped  nodes ) ; 

%%  End  Assembly 

fprintf ( fid2 , ' *End  AssemblyVn ' ) ;  %ends  the  assembly  portion  of  it 
fprintf ( f id2 , ' **\n ' ) ; 

%%  Materials 
%Standard  Heading 
fprintf (fid2,  ' **  MATERIALS \n ' )  ; 
fprintf ( f id2 , ' **\n ' ) ; 

%  Determine  and  Assign  Matl  Prop  for  the  Veins  (B32  elements) 
for  x=l : size (angle, 1 ) 

theta  =  angle (x, 2 ) +theta  original ( 1 ,:) ; 

%set  up  material  properties 

fprintf (fid2,  ' *Material,  name=Material-%-2 .  Of \n ' , x) ; 

fprintf ( f id2 , ' *Density\n ' ) ; 

fprintf (fid2, ; %#6 . 6f\n ' , Density (x) ) ; 

fprintf (fid2, ' *Elastic,  type=ENGINEERING  CONSTANTS \n ' ) ; 

%  Use  the  material  propertie  saved  from  Section  section 
%  determine  which  properties  to  use 
ctr=max (Matl_Ex (x,  : ) )  ; 
if  Matl_Ex(x,l)  ==  ctr; 
ctr  =  1; 

elseif  Matl_Ex(x,2)  ==  ctr; 
ctr=2 ; 

elseif  Matl_Ex(x,3)  ==  ctr; 
ctr=3 ; 

end 

fprintf (fid2, ' %#12 . 6g,  %#12.6g,  %#12.6g,  %#12.6g,  %#12.6g,  %#12.6g, 
%#12 . 6g,  %#12 . 6g, \n ' . . . 

,  Matl_Ex (x,  ctr ) , Matl_Ey (x, ctr ) , Matl_Ey (x, ctr ) ,0.1, 0.1, 0.1,  .  .  . 
Matl_Gxy (x, ctr) , Matl_Gxz (x, ctr) ) ; 
fprintf (fid2, ;%#12.6g,\n', Matl_Gyz (x, ctr ) ) ; 

end 


%%  Step 

%Standard  Heading 

fprintf ( fid2 , '**  - 

- \n') ; 

fprintf (fid2,  ' **  \n'); 

fprintf (fid2, ' **  STEP:  Step-l\n'); 

fprintf ( f id2 , ' **\n ' ) ; 

%since  there  is  only  one  step  and  we  always  want  to  find  freqency, 
it'll  look  like  this 

fprintf ( fid2 ,' *Step,  name=Step-l,  perturbation\n ' ) ; 
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fprintf ( f id2 , ' Find  Freq.  \n'); 

fprintf ( f id2  ,  ' *Frequency,  eigensolver=Lanczos ,  acoustic  coupling=on, 

normalization=displacement\n ' )  ; 

fprintf (fid2,  ' %6. Of ,  ,  ,  ,  , \n ' , num_modes ) ; 

%%  Boundry  Conditions 
%  Standard  Heading 
fprintf ( f id2 , ' **  \n '  )  ; 

fprintf (fid2, ' **  BOUNDARY  CONDITIONS\n ' ) ; 
fprintf (fid2,  '  **  \n  '  )  ; 

%create  a  boundry  condition  (in  this  case  it  will  be  fixed) 
fprintf (fid2 ,' **  Name:  BC-1  Type:  Displacement/Rotation\n ' ) ;  %since 
only  one,  can  fix  6  dof,  says  type  of  boundry  cond 


fprintf 

( f id2 , 

' *Boundary\n ' ) ; %tells 

abaqus  there  is 

a  boundry 

conditon 

fprintf 

dir 

( f id2 , 

'  PickedSet%-2 . Of , 

1, 

l\n  ' 

, (set  ctr+3) 

) ; %hold 

in 

x  or  1 

fprintf 

dir 

( f id2 , 

'  PickedSet%-2 . Of , 

2, 

2\n ' 

, (set  ctr+3) 

) ; %hold 

in 

y  or  2 

fprintf 

dir 

( f id2 , 

'  PickedSet%-2 . Of , 

3, 

3\n ' 

, (set  ctr+3) 

) ; %hold 

in 

z  or  3 

fprintf 
1  axis 

( f id2 , 

'  PickedSet%-2 . Of , 

4, 

4 \n  ' 

, (set  ctr+3) 

) ; %hold 

in 

rot  about 

fprintf 
2  axis 

( f id2 , 

'  PickedSet%-2 . Of , 

5, 

5\n ' 

, (set  ctr+3) 

) ; %hold 

in 

rot  about 

fprintf 
3  axis 

( f id2 , 

'  PickedSet%-2 . Of , 

6, 

6\n ' 

, (set  ctr+3) 

) ; %hold 

in 

rot  about 

%%  Outputs 
%Standard  Heading 
fprintf ( f id2  ,  ' **\n ' ) ; 

fprintf (fid2, ' **  OUTPUT  REQUESTS \n ' ) ; 
fprintf ( f id2 , ' **\n ' ) ; 

%Frequency  outputs  (in  this  case  it  does  not  change) 
fprintf (fid2, ' *Restart,  write,  frequency=0 \n ' ) ; 
fprintf ( f id2 , ' **  \n '  )  ; 

fprintf (fid2, ' **  FIELD  OUTPUT:  F-Output-l\n ' ) ; 
fprintf ( f id2 , ' **  \n '  )  ; 

fprintf ( fid2 ,' *Output ,  field,  variable=PRESELECT\n ' ) ; 

%%  End  Step 

fprintf ( fid2 ,' *End  Step\n'); 

%%  Close  the  File 
fclose ( f id2 ) ; 

%%  Run 

%  This  line  runs  the  completed  file 

eval ( [ ' dos ( ' ' abq6102  j ob=Beam_Test_out . inp  interactive'')']) 

%%  Open  . dat  File 

%Gather  the  results  from  the  .dat  file  and  store  in  a  matrix 
fid3=fopen ( [File  '  . dat ' ]  ,  ' r+ ' )  ;  %  Open  file  for  reading/writing 
idx  =  0 ; 

while  ~feof(fid3) 
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OUTPUT  '  )  ) 


idx  =  idx  +  1; 
tline=fgetl (fid3) ; 

if  isempty (strfind (tline, 'EIGENVALUE 
==  0  %  Scan  till  *Nodes 
for  i=l : 5 

tline=f getl ( f id3 ) ;  %  Read  5  lines  of  blank  stuff 

end 

Eigen  value_string=fgetl ( f id3 ) ; %calculate  the  new  value  of 
maximum  nodes 

make  text  =  str2num  (Eigenvalue  string);  %converts  text  into 
numbers  MATLAB  can  read 

Eigen  value  =  make  text (1,4);  %  gets  the  Frequency  from  the 
string  of  numbers 
end 

end 

fclose ( f id3 ) ;  %  Close  the  File 
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Appendix  D.  Beam  .inp  File 


This  is  the  original  .inp  file  that  was  solved  for  the  wire  beam  element. 

*Heading 

**  Job  name:  BeamTest  Model  name:  Model-1 
**  Generated  by:  Abaqus/CAE  6.10-2 
*Preprint,  echo=NO,  model=NO,  history=NO,  contact=NO 
** 

**  PARTS 
** 

*Part,  name=Part-l 
*End  Part 

** 

**  ASSEMBLY 
** 

*Assembly,  name=Assembly 
** 

*lnstance,  name=Part-l-l,  part=Part-l 
*Node 


1, 

0., 

0., 

0. 

2, 

4., 

0., 

0. 

3, 

8., 

0., 

0. 

4, 

12., 

o., 

0. 

5, 

16., 

o., 

0. 

6, 

20., 

o., 

0. 

7, 

24., 

o., 

0. 

8, 

28., 

o., 

0. 

9, 

32., 

o., 

0. 

10, 

36., 

0., 

0. 

11, 

40., 

0., 

0. 

12, 

2., 

0., 

0. 

13, 

6., 

0., 

0. 

14, 

10., 

0., 

0. 

15, 

14., 

0., 

0. 

16, 

18., 

0., 

0. 

17, 

22., 

0., 

0. 

18, 

26., 

0., 

0. 

19, 

30., 

0., 

0. 

20, 

34., 

0., 

0. 

21, 

38., 

0., 

0. 

*Element,  type=B32 

1,  1,12,  2 

2,  2,  13,  3 

3,  3,  14,  4 

4,  4,  15,  5 

5,  5,  16,  6 

6,  6,  17,  7 

7,  7,  18,  8 

8,  8,  19,  9 

9,  9,  20,  10 

10,  10,21,  11 

*Nset,  nset=_PickedSet2,  internal,  generate 
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1,  21,  1 

*Elset,  elset=_PickedSet2,  internal,  generate 

1,  10,  1 

*Nset,  nset=_PickedSet4,  internal,  generate 

1,  21,  1  ' 

*Elset,  elset=_PickedSet4,  internal,  generate 

1,  10,  1 

*Orientation,  name=Ori-l 

1.. 0,  0.,  0.,  1.,  0. 

1,0. 

**  Section:  Section- 1  Profile:  Profile-1 

*Beam  Section,  elset=_PickedSet2,  material=Material-l,  temperature=GRADIENTS,  section=I 
0.5,  1.,  2.5,  2.5,  0.33,  0.33,  1. 

0.,0.,-l. 

*End  Instance 

** 

*Nset,  nset=_PickedSet4,  internal,  instance=Part-l-l 
1, 

*End  Assembly 
** 

**  MATERIALS 
** 

*Material,  name=Material-l 
*Density 

1785., 

*Elastic 
4.15e+09,  0.28 

** _ 

** 

**  STEP:  Step-1 
** 

*Step,  name=Step-l,  perturbation 

*Frequency,  eigensolver=Lanczos,  acoustic  coupling=on,  normalization=displacement 

10. .  . . , 

H=H= 

**  BOUNDARY  CONDITIONS 

** 

**  Name:  BC-1  Type:  Displacement/Rotation 

*Boundary 

_PickedSet4,  1,  1 

_PickedSet4,  2,  2 

_PickedSet4,  3,  3 

_PickedSet4,  4,  4 

_PickedSet4,  5,  5 

_PickedSet4,  6,  6 

H=H= 

**  OUTPUT  REQUESTS 
** 

*Restart,  write,  frequency-0 

** 

**  FIELD  OUTPUT:  F-Output-1 

H=H= 
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Appendix  E.  Fiber  Angle  Measurement  MATLAB  Code 

This  is  the  MATLAB  code  developed  to  measure  the  angle  of  the  fibers  in  the 
images  taken  by  the  optical  microscope  of  the  carbon  composite. 

%Find  the  angle  of  the  fiber  orientation  vs  the  edge 

%  click  along  the  edge  first,  then  click  5  fibers  (10  clicks)  on  random 
%  fibers  to  find  the  average  fiber  orientation  of  the  specimen 
clc, clear, close  all 

%%  First,  let's  read  all  the  data-sets  from  the  damO.unv  file 
[FileName, PathName]  =  uigetf ile ( ' * . tif * '  ,  ' 0  ' )  ; 

File  =  FileName (1 : findstr (FileName,  '.') -1)  ; 

Ext  =  FileName ( (findstr (FileName, '.') +1) : length (FileName) ) ; 

%%  Read  in  Moth  image 
I  =  imread ([ PathName  FileName]); 

%%  Plot  the  picture  in  a  figure 
imshow ( I ) ; 

%%  Pick  3  Points  on  the  image 
%  location  is  in  pixels 

warndlg ( ' Pick  the  first  set  of  six  Fibers !!','!!  Warning  ! ! ' ) 

[x, y] =ginput (12) ; 

%%  Save  .mat  file  with  image  filename 
save ( [File  ' .mat ' ] ) 

%first  fiber 

Theta_Fiberl=abs (atand ((y(2)-y(l))/(x(2)-x(l))));%pick  one  orientation 
first 

Theta_Fiber2=abs (atand ( (y ( 3 ) -y ( 4 ) ) / (x ( 3 ) -x ( 4 )  )  )  )  ;  %  pick  one  fiber 
Theta_Fiber3=abs (atand ( (y ( 5 ) -y ( 6 ) ) / (x ( 5 ) -x ( 6 ) ) ) ) ; %  pick  second  fiber 
Theta_Fiber4=abs (atand ( (y ( 7 ) -y ( 8 ) ) / (x ( 7 ) -x ( 8 ) ) ) ) ; %  pick  third  fiber 
Theta_Fiber5=abs (atand ( (y ( 9 ) -y ( 10 ) ) / (x ( 9 ) -x ( 10 ) ) ) ) ; %  pick  fourth  fiber 
Theta_Fiber6=abs (atand ( (y (11) -y (12) ) / (x (11) -x (12) ) ) ) ; %  pick  fifth  fiber 
Theta  Fiber  avgl  =  (Theta  Fiberl  +  Theta  Fiber2  +  Theta  Fiber3  + 

Theta  Fiber4 . . . 

+  Theta  Fiber5 ) /5 ; %f ind  average 

warndlg ( ' Pick  the  second  set  of  six  Fibers !!','!!  Warning  ! ! ' ) 

[x, y] =ginput (12) ; 

%second  fiber 

Theta_Fiberl=abs (atand ((y(2)-y(l))/(x(2)-x(l))));%pick  one  orientation 
first 

Theta_Fiber2=abs (atand ( (y ( 3 ) -y ( 4 ) ) / (x ( 3 ) -x ( 4 ) ) ) )  ;  %  pick  one  fiber 
Theta_Fiber3=abs (atand ( (y ( 5 ) -y ( 6 ) ) / (x ( 5 ) -x ( 6 ) ) ) ) ; %  pick  second  fiber 
Theta_Fiber4=abs (atand ( (y ( 7 ) -y ( 8 ) ) / (x ( 7 ) -x ( 8 ) ) ) ) ; %  pick  third  fiber 
Theta_Fiber5=abs (atand ( (y ( 9 ) -y ( 10 ) ) / (x ( 9 ) -x ( 10 ) ) ) ) ; %  pick  fourth  fiber 
Theta_Fiber6=abs (atand ( (y (11) -y (12) ) / (x (11) -x (12) ) ) ) ; %  pick  fifth  fiber 
Theta  Fiber  avg2  =  (Theta  Fiberl  +  Theta  Fiber2  +  Theta  Fiber3  + 

Theta  Fiber4 . . . 

+  Theta  Fiber5 ) /5 ; %f ind  average 
Theta  =  abs (Theta  Fiber  avg2-Theta  Fiber  avgl) 
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Appendix  F.  Beam  Monte  Carlo  Code 


This  developed  MATLAB  code  was  used  to  submit  multiple  .inp  fdes  used  in  the 

Monte  Carlo  Solution  to  MATLAB.  Data  was  then  stored  in  a  matrix  and  plotted  when 

the  solution  was  finished  running. 

clc, clear  all;  close  all; 

%%  Define  Number  of  Cases  to  Run 
NumCases=7500; 

%%  Definite  Variation  in  Laser  to  Fiber  Orientation 
pm=5;  %This  represents  +/- 
nstd  alpha=l; 
avg  alpha=2.5; 

alpha=avg  alpha+pm*rand (NumCases ,  1 ) -pm/ 2  ; 
alpha= [alpha  alpha  alpha]; 

%%  Define  Variation  in  Lamina  Layup  Orientation 
pm=10;  %This  represents  +  /- 
nstd_theta=l ; 
avg_theta=2 . 5 ; 

theta  original  =  avg_theta+pm*rand (NumCases , 3 ) -pm/2 ; 
theta_original ( : , 2 ) =90-theta_original ( : , 2 ) ; 

%%  Define  Variation  in  Laminate  Thickness 
pm=10 ; 
nstd_t=l ; 

avg_thickness=145; 

h  =  avg  thickness+pm*rand (NumCases , 1 ) -pm/2 ; 
h  =  h*lE-6;  %convert  from  integers  to  microns 
%%  Define  Parameter  Estimates 

[ahat, bhat, ACI , BCI ]  =  unifit(theta  original, . 01) ; 

%%  Run  the  Analysis 
for  t=l: NumCases 

[Eigenvalue  (t,  :  )  ]  =Vary_Properties_Beam  (alpha  (t)  ,  theta_original  (t,  :  )  ,  h  ( 

t)  , t)  ; 

end 


%%  Plot  the  Results  Histogram  -  Modal  Frequency 
close  all 

%  Open  a  Figure  the  size  of  the  screen 
scrsz  =  get (0,  ' ScreenSize ' ) ; 

figure (' Position [1  scrsz (2)  scrsz (3)  scrsz (4) ] ) 

nbins=100 ; 
subplot  (331) 
hist (alpha, nbins) 
grid  on 

xlabel ([' Alpha  =  Mean  '  num2str (avg_alpha)  '  deg  |  std  = 
num2str(nstd  alpha)]) 
ylabel ( ' Bins ' ) 

subplot  ( 332 ) 
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hist (theta  original (:, 1 ), nbins ) 
grid  on 

xlabel ([' Theta  =  Mean  '  num2str (avg_theta)  '  deg  |  std  =  ' 

num2str (nstd_theta) ] ) 
ylabel ( ' Bins ' ) 

subplot ( 333 ) 
hist(h(:,l)*lE6, nbins) 
grid  on 

xlabel ([' Thickness  =  Mean  '  num2str(avg  thickness)  '  um  |  std  =  ' 

num2str (nstd_t) ] ) 
ylabel ( ' Bins ' ) 

subplot ( 3 , 3 , 4 : 6 ) 

hist (Eigen  value, nbins) 

grid  on 

xlabel (' Frequency  (Hz) ') 
ylabel ( ' Bins ' ) 

title (['Mean  freq  =  '  num2str (mean (Eigen_value) , ' %f ' )  '  (Hz)']) 

subplot ( 3 , 3 , 7 : 9 ) 

plot  (1  :NumCases,  Eigenvalue,  '  .  ,  'MarkerSize  ,  8) 
grid  on 

xlabel ( ' Sample ' ) 
ylabel (' Frequency  (Hz) ') 
hold  on 

plot ( [0  NumCases ] , mean (Eigen  value) * [ 1  1] , ' — g ' , ' LineWidth ' , 2 ) 
legend ( ' Samples ' , ' Mean ' , 'Location ; , ! SouthEast ' ) 

title (['Mean  freq  =  '  num2str (mean (Eigen  value), ' %f ' )  '  (Hz),  Max  freq 

=  '  num2str (max (Eigen  value),  ' % f ' )  '  (Hz),  Min  freq  = 

num2str (min (Eigenvalue) <  ' %f  '  )  '  (Hz)']) 

hold  off 
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Appendix  G.  Wing  Develop  .inp  File  Code 


This  is  the  code  that  was  used  to  develop  the  FEA  model  of  the  Engineered  wing. 


function [Eigen  valuel, Eigen  value2 ] =Final  Wing (alpha, theta_original,h, t 

) 

%function-Outputs  =  First  and  Second  Modal  Frequency 

Inputs  =  laser  cut  angle,  ply  orientation,  thickness,  run  number 
%File  names 

Input  File  =  'Wing  Shoulder  Eng  Doe. inp';  %  Name  of  the  file  created 
with  extension 

Node  File  =  'WING  FIT  SHOULDER  DD.inp';  %  Name  of  file  where  nodes  are 
located  with  extension 

%Scale  Nodes 

%Note,  make  root  of  wing  =  0,0. 

x  offset=0;  %Move  nodes  in  x  dir 

y_offset=0 . 592322826000000;  %Move  nodes  in  y  dir 

wing^length  =  0.05;  %Desired  length  of  wing  in  meters 

%Scale  of  the  wing 

%Begining  and  ending  width  of  the  vein,  linearly  tapering  in  meters  x 
10^5 

Costa  =  [75  20];  %Costa  Vein 
Radius  =  [80  60];  %  Radius  Veins 
Archulus  =  [60  35];%  Archulus 
RMCA  =  [55  15];  %  Rest  of  the  veins 

StringersWidth= [Costa;  %  Stringer  1  -  Leading  Edge  Vein 

Costa;  %  Stringer  2 

Costa;  %  Stringer  3 

Radius;  %  Stringer  4  -  Medial  Vein 

Archulus;  %  Stringer  5  -  Arculus  Vein 

RMCA;  %  Stringer  6 

RMCA;  %  Stringer  7 

RMCA;  %  Stringer  8 

RMCA;  %  Stringer  9 

RMCA;  %  Stringer  10 

RMCA;  %  Stringer  11 

RMCA; ] *lE-5;  %  Stringer  12 

%Material  Properties 
%Carbon  fiber  -  IDEAL  VALUES 

%  alpha  =  0;  %Angle  misallignment  placed  in  laser 
%  theta_original=  [0  90  0];  %Fiber  orientation 
%  h  =  150e-6;  %Thickness 

%Material  Properties 
%Carbon  fiber 

hi  =  h;  %Thickness  part  2,  must  be  same  as  h  (used  to  reset  h  in  some 
calc) 

den  cf  =  1790*0.88;%  Density  of  the  carbon-fiber  composite 
%Membrane 

shell  thick  =  20.0e-6;  %Thickness  of  the  membrane  (um) 
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den  sh  =  1420*0.88;  %Density  of  the  Mylar  membrane  (kg/m^3) 
shell  E  =  2.5E9;  %Modulus  of  Membrane  (GPa) 
shell  v  =  0.3;  %Poisson's  Ratio  of  Membrane 

%Nodes 

Base  Node  Loc=  [2,48,51,280,283,942,1269];  %Node  number  where  you  want 
to  apply  the  Boundry  Condition. 

Pt  Node  Loc  =  1108;  %  Node  number  where  Pt  load  is  applied 
%Frequency 

num  modes  =  9;  %The  number  of  modes  you  wish  to  find 
%Point  Load 

Pt  Load  =  -0.001;  %The  force  you  want  to  put  on  the  tip  of  the  wing 
%%  Open  the  File 

fid2=fopen (Input_File, 'w');  %  Open  file  for  writing,  new  file 
%%  Heading 

%  Standard  Abaqus  Heading 
fprintf ( f id2  ,  ' * Heading \n ' ) ; 

fprintf ( f id2 , '**  Job  name:  BEAM  MODAL  BASELINE  Model  name:  %s\n',date) 
%%  Contact 

%  Contact  within  model 

fprintf (fid2,  ' *Preprint,  echo=NO,  model=NO,  history=NO,  contact=NO\n ' ) 
fprintf ( f id2 , ' **\n ' ) ; 

%%  Parts 

%  Parts  (doesn't  change  for  this  case) 
fprintf (fid2,  ' **  PARTS\n'); 

fprintf (fid2, ' *Part,  name=WING  FIT  SHOULDER  EE  %s\n',date);  %  name 

called  up  in  instance 

fprintf ( f id2 , ' **\n ' ) ; 

fprintf ( fid2 ,' *End  Partin'); 

fprintf ( f id2 , ' **\n ' ) ; 

fprintf (fid2, '**'); 

%%  Assembly 

%insert  the  Assembly  section  (doesn't  change  for  this  case) 
fprintf (fid2, ' **  ASSEMBLY\n ' ) ; 
fprintf (fid2, ' **\n' ) ; 

fprintf ( fid2 , ' *Assembly,  name=Assembly\n ' ) ; 
fprintf ( f id2 , ' **\n ' ) ; 

%%  Instance 

%  Insert  the  Instance  section 

fprintf ( fid2  ,  ' *Instance,  name=Part-l-l , 

part=WING  FIT  SHOULDER  EE  %s\n '  ,  date) ; %  name  called  from  parts 

%%  Insert  the  nodes  from  the  veins 
%Standard  Heading 

fprintf ( fid2 ,' *Node ') ;  %heading  for  the  nodes 
fprintf ( f id2 , ' \n ' ) ; 
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NodesMax  =  0;  %Placeholder  for  number  of  nodes  counter  (cacl  later) 


%  Take  the  nodes  out  of  another  inp  file. 

fidl=fopen (Node  File,  r+ ' )  ;  %  Open  file  for  reading/writing 
idx  =  0 ; 

while  ~feof(fidl) 
idx  =  idx  +  1; 
tline=fgetl (fidl) ; 

if  isempty (strfind (tline, ' *Node ' ) )  ==  0  %  Scan  till  *Nodes 

for  i=l : 2000000000  %  Scan  till  *  Found 
tline=fgetl (fidl) ; 

if  isempty (strfind (tline, '*') )  ==  0 

NodesMax=i-l ; %calculate  the  new  value  of  maximum  nodes 
break 

else 

nodes (i,  : ) =sscanf (tline,  %d,%f,%f,%f')  '  ; 

end 

end 

break 

end 

end 

fclose (fidl ) ; 


%  for  loop  to  print  nodes  into  the  .inp  file 
%  nodes  will  be  scaled  at  this  time 

node_scale  =  wing_length/ (max (nodes ( : , 2 ) ) ) ;  %sets  scale  =  to 
legnth 

fprintf ( f id2  ,  ' %d,  %2.10f,  %2.10f,  %2 . 10f\n ',  [nodes (:, 1 ),..  . 

(nodes ( : , 2 ) - 


x_offset) *node_scale,  .  .  . 
y_offset) *node_scale,  .  .  . 


(nodes ( : , 3) - 
nodes ( : , 4 ) ] ' ) ; 


desired 


%%  Elements 
%  For  beam  elements 

fprintf (fid2,  ' *Element,  type=B32 ' ) ;  %heading  for  the  elements 
fprintf ( f id2 , ' \n ' ) ; 

%Search  Input  File  for  *Elem  Locations 

fidl=fopen (Node  File,  ; r+ ' )  ;  %  Open  file  for  reading/writing 
idx  =  0 ; 

while  ~feof(fidl) 
idx  =  idx  +  1; 
tline=fgetl (fidl) ; 

if  isempty (strfind (tline, ' *Element,  type=B32 ' ) )  ==  0  %  ' ' 
for  i=l: NodesMax  %  Scan  till  *  Found 

tline=fgetl (fidl) ; 

if  isempty (strfind (tline,  '*') )  ==  0 
break 

else 

elems . B32 (i,  :) =sscanf (tline,  '%f,%f,%f,%f,%f,%f,%f,%f,%f') 
end 

end 
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break 


end 

end 

fclose ( f idl ) ; 

%  Write  nodes  into  the  program  file 

fprintf ( f id2 ,  ' %d, %6 . Of , %6 . Of , %6 . Of \n ' ,  [elems.B32]  '  )  ; 

%  Heading  for  quad  elements  membrane  section 

fprintf (fid2, ' *Element,  type=S8R');  %heading  for  the  elements 
fprintf ( f id2 , ' \n ' ) ; 


%Search  Input  File  for  *Elem  Locations 
fidl=fopen (Node  File, 'r+'); 
idx  =  0 ; 

while  ~feof(fidl) 
idx  =  idx  +  1; 
tline=fgetl (fidl) ; 

if  isempty (strfind (tline, ' *Element,  type=S8R'))  ==  0  %  Scan  till  '' 
for  i=l:NodesMax 

tline=fgetl (fidl) ; 

if  isempty (strfind (tline, '*') )  ==  0 
break 

else 


elems.S8R(i, :) =sscanf (tline, 
end 


»  O,  -p  0,4=  0,4=  0,4=  O,  4=  0,4=  0 .4=  0,4=  O.  r  |  \  |  . 

oi  ,  ol,  o  L  f  ol  ,  o  L  f  ol  ,  o  L  f  O  L  f  oJ_  ;  , 


end 

break 


end 


end 

fclose (fidl ) ; 


%Write  the  Elements  to  the  File 

fprintf (fid2,  ' %d,  %6 . Of , %6 . Of , %6 . Of , %6 . Of , %6 . Of , %6 . Of , %6 . Of , %6 . Of \n ’ ,  [el 

ems . S8R (:,:)]'); 


%  For  tri  elements  membrane  section 

fprintf (fid2, ' *Element,  type=STRI 65 ' ) ;  %heading  for  the  elements 
fprintf ( f id2 , ' \n ' ) ; 


%Search  Input  File  for  *Elem  Locations 

fidl=fopen (Node  File,  '=+');  %  Open  file  for  reading/writing 
idx  =  0 ; 

while  ~feof(fidl) 
idx  =  idx  +  1; 
tline=fgetl (fidl) ; 

if  isempty (strfind (tline, ' *Element,  type=STRI 65 ' ) )  ==  0  %  Scan  till 
for  i=l:NodesMax 

tline=fgetl (fidl) ; 

if  isempty (strfind (tline,  '*') )  ==  0 
break 

else 


elems . S TRI 65 (i, :)=sscanf (tline, 
end 


»  9-  f  9--F  9--F  9--F  9--F  S-  -F  9--F  9-  -F  9-f  M  '  • 

O  J -  f  ol  ,  ol  ,  O-Lj  ol  ,  ol,  ol  ,  o  )  f 


end 
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break 


end 

end 

fclose ( f idl ) ; 

%  Print  Tri  Elements  to  the  Program 

fprintf(fid2,  '%d,%6.0f,%6.0f,%6.0f,%6.0f,%6.0f,%6.0f\n', elems . STRI 65 ' ) ; 

%%  Read  Stringer  Node  Sets 

%Search  Input  File  for  Stringer  Nsets 

fidl=fopen (Node  File, ' r+ ' ) ;  %  Open  file  for  reading/writing 
for  j  =1 : 12  % (Number  of  Stringers) 
idx  =  0; 

Stringers  ( j )  . nset=  [  ] ; 
while  ~feof(fidl) 
idx  =  idx  +  1; 
tline=fgetl (fidl) ; 

if  isempty (strfind (tline, ['*Nset,  nset=  Stringer-'  num2str(j) 

' ,  internal ']))==  0  %  ' ' 

for  i=l:NodesMax  %  Scan  till  *  Found 

tline=fgetl (fidl) ; 

if  isempty (strfind (tline, '*') )  ==  0 
break 

else 

Stringers (j) . nset= [ Stringers ( j ) .nset 
sscanf (tline, ' %f , '  ) '  ]  ; 

end 

end 

break 

end 

end 

end 

fclose (fidl ) ; 

%%  Read  Stringer  Element  Sets 
%Search  Input  File  for  Stringer  Elsets 

fidl=fopen (Node  File,  ' r+  '  )  ;  %  Open  file  for  reading/writing 
for  j  =1 : 12  % (Number  of  Stringers) 
idx  =  0; 

Stringers  ( j )  . elset= [ ] ; 
while  ~feof(fidl) 
idx  =  idx  +  1; 
tline=fgetl (fidl) ; 

if  isempty (strfind (tline, [' *Elset,  elset=_Stringer- '  num2str(j) 
',  internal,  generate']))  ==  0  %  '' 

for  i=l:NodesMax  %  Scan  till  *  Found 

tline=fgetl (fidl) ; 

if  isempty (strfind (tline, '*') )  ==  0 
break 

else 

Stringers (j) . elset= [ Stringers ( j ) .elset 
sscanf (tline, ' %f , ' ) ' ] ; 

end 

end 

break 

end 

end 
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end 

fclose ( f idl ) ; 

%%  Process  Stringer  Element  Sets 
for  j  =1 : 12 

Stringers(j)  . vect= [ Stringers  ( j )  .elset(l)  : Stringers ( j )  .elset(3)  :Stringer 
s ( j ) .elset (2) ] ' ;  %Vector  of  Elements 

%  Find  Index  of  Stringer  Elements  within  the  subset  of  all  the 
Model  Elements 

[c,  ia,  ib] =intersect (elems .B32 (:, 1) , Stringers (j ) .vect) ;  %  c  = 
common  values,  ia  =  index  values  in  A,  ib  =  index  values  in  b 

%  Assign  Stringer  Nodes  to  Elements 
Stringers ( j )  . elems=elems . B32 ( ia, 1 : end)  ; 

%  Assign  Midpoint  Node  Location  to  Elems 
Stringers (j) . elems= [ Stringers ( j ) .elems 
nodes (Stringers ( j ) . elems ( : , 3 ) , 2 : end) ] ; 

%  Sort  Elems  based  on  Midpoint  Nodes  X-Location 
if  j  ~=  5 

[ Stringers ( j ). elems ] =sortrows ( Stringers ( j ). elems , 5 ) ;  %  Sort 
veins  based  on  X-Location 
else 

[ Stringers ( j ). elems ] =sortrows ( Stringers ( j ). elems , 6 ) ;  %  Sort 
Arculus  based  on  Y-Location 
end 

end 

for  j  =1 : 12 

Stringers ( j ) . widths=linspace ( StringersWidth ( j , 1) , StringersWidth ( j , 2) , si 
ze (Stringers ( j ) .elems, 1) ) 

[Stringers ( j )  . elems ]  =  [Stringers ( j )  . elems  Stringers ( j )  . widths ]  ;  % 
Assign  widths  to  last  column 
end 


%%  Rearrange  Structured  Stringers  Variable  into  new  Unstructured 
Variable 

elems . Stringers= [] ;  %  [Elem  #,  StartNode, MidNode, EndNode, Width,  (add 
angle  later  in  code) ] 
for  j  =1 : 12 

elems . Stringers= [elems . Stringers ;  Stringers ( j ) . elems ( : , [1:4,8])]; 

end 

elems . Stringers=sortrows (elems . Stringers , 1 ) ;  %  Sort  Elems  by  Elem# 

%%  Calculate  Angle  between  the  Stringer  B32  Elements 
%  needed  to  determine  material  properties, 
for  m  =  1 : (length (elems . B32 ( : , 1) ) ) 

%determine  the  first  node  that  is  part  of  the  element  i.e  x 
coordinate 

node  str  =  elems . B32 (m, 2 )  ; 
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%determine  the  second  (mid  point)  node  of  the  element  i.e  y 
coodiinate 

node  end  =  elems . B32 (m, 3 ) ; 

%location  of  first  node 
loc_str  = 

[ (nodes (node_str , 2 ) ) , (nodes (node_str , 3 ) ) , (nodes (node^str, 4 ) ) ] ; 

%  location  of  the  end  node 
loc_end  = 

[ (nodes (node_end, 2 ) ) , (nodes (node_end, 3 ) ) , (nodes (node_end, 4 ) ) ] ; 

%  determine  the  slope  of  the  line 

slope_2  =  ( loc  str ( 1 , 2 )  -  loc_end ( 1 , 2 ) ) / (loc_str (1 , 1 )  - 
loc_end (1,1) ) ; 

%determine  the  angle  with  horizontal  [inf, 0,0]  vector  from  origin 
theta_elem  =  atand ( slope_2 ) ; 

%store  it! 

angle (m, 1 ) =m;  %stores  the  row  number  in  angle  in  1st  column 
angle (m, 2) =  elems . B32 (m, 1 ) ;  %the  2nd  column  hold  element  number 
angle (m, 3 ) =theta_elem;  %  3rd  column  holds  the  angle 

end 

elems . Stringers= [elems . Stringers  angle(:,3)];  %  Append  Angle  Data 

%%  Write  Nset  and  Elsets 
%  For  the  beam 

set  ctr=0;  %counts  the  number  of  Node  and  Element  Sets 
%  Create  a  new  section  for  each  Vein  element 
for  x=l : length (angle) 

%create  heading  and  pick  a  Nset 
fprintf (fid2, ' *Nset,  nset=  PickedSet%-2 . Of , 
internal \n  , angle (x, 1 ) ) ; 

nset  nodes  =  sort (elems . B32 (x, 2 : 4 ) ) ; %sort  out  the  nodes  in 
asscending  order 

%print  nodes  in  that  Nset  i.e.  for  each  element 
fprintf (fid2, '%6. Of,  %6.0f, 

%6 . 0f\n ' , nset_nodes (1) , nset_nodes (2) , nset_nodes (3) ) ; 

%  for  Elsets,  create  elset  same  number  as  corresponding  Nset 
fprintf (fid2, ' *Elset,  elset=  PickedSet%-2 . Of , 
internal \n ' , angle (x, 1 ) ) ; 

fprintf (fid2, ' %6 . Of \n ' , angle (x, 2 ) ) ; 

set_ctr=set_ctr+l ;  %add  to  Set  counter 

end 

%  For  the  Membrane 
%  Quad  S8R  elements 

%  Sort  the  nodes  in  order,  using  unique  to  get  rid  of  duplicates 
s8r_node=sort (unique (elems . S8R) ) ; 

fprintf ( fid2  ,' *Nset,  nset=  PickedSet%-2 . Of ,  internal\n '  ,  (set  ctr+1)  )  ; 
fprintf (fid2 ,  ' %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d, 
%d,  %d\n ' , s8r_node ( :  ,  1 ) )  ; 
fprintf ( f id2 , ' \n ' ) ; 

%  sort  the  elements  in  order 
s8r_elem=sort (elems . S8R ( :  ,  1 ) )  ; 

%  for  Elsets,  create  elset  same  number  as  corresponding  Nset 
fprintf (fid2, ' *Elset,  elset=  PickedSet%-2 . Of ,  internal\n ' , (set  ctr+1)) 
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fprintf (f id2 ,  1 %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d, 
%d,  %d\n ' , s8r_elem ( : , 1) )  ; 
fprintf ( f id2 , ' \n ' ) ; 

%  Tri  STRI65  elements 

%  Sort  the  nodes  in  order,  using  unique  to  get  rid  of  duplicates 
tri_node=sort (elems . STRI 65 )  ; 

fprintf (fid2,  ' *Nset,  nset=  PickedSet%-2 . Of ,  internal\n ' ,  (set  ctr+2)); 
fprintf ( fid2 ,  ' %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d, 

%d,  %d\n ' , tri_node ( :  ,  1 )  )  ; 
fprintf ( f id2 , ' \n  '  )  ; 

%  sort  the  elements  in  order 
tri_elem=sort ( (elems . STRI 65 ( :  ,  1 ) ) )  ; 

%  for  Elsets,  create  elset  same  number  as  corresponding  Nset 
fprintf (fid2, ' *Elset,  elset=  PickedSet%-2 . Of ,  internal\n ' , (set  ctr+2)); 
fprintf (fid2 ,  ' %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d,  %d, 

%d,  %d\n' , tri_elem ( : , 1) ) ; 
fprintf ( f id2 , ' \n  '  )  ; 

%%  Beam  Sections 

%  The  sections  will  start  with  the  veins  and  begin  with  Section  1  (B32) 

for  x=l : length (angle)  %  new  section  for  each  element 

fprintf ( fid2 ,'* *  Section:  Section-%-2 . Of  Profile:  Profile-%- 

2 . Of \n ' . . . 

,x,x);  %  create  section/profile 
fprintf (fid2, . . . 

' *Beam  Section,  elset=  PickedSet%-2 . Of ,  material=Material-%- 
2. Of,  temperature=GRADIENTS ,  section=I\n ' . . . 

,x,x);  %  match  the  element  with  it's  own  elset  and  material 

b (x) =elems . Stringers (x, 5) ; %artificial  width  of  the  beam  element 

%  Generate  the  I  beam  sections 
%  Generate  the  I  beam  sections 

theta  =  angle (x, 3 ) +theta  original ( 1 ,:) ;  %determine  theta 
%For  the  top  ply  run  a  function  to  determine  Matl  Prop 
[ ~ , Ex, Ey, ~ , ~ , Vxy, Vyx, Gxy, Gxz , Gyz  ,  ~ ,  ~ ,  ~ ,  ~ ]  =  .  .  . 

LAMINATE_NASA_HALPIN_HYBRID_BASELINE_FUNC (4  0E- 
3, b (x) , h, alpha, theta (1,1)  ,  0 )  ; 

%  Capture  the  Mat  Properties  so  we  don't  have  to  run  this  again 
Matl  Ex(x,l)  =  Ex;  Matl  Ey(x,l)  =  Ey;  Matl  Vxy(x,l)  =  Vxy; 

Matl  Vyx(x,l)  =  Vyx; 

Matl_Gxy (x, 1 )  =  Gxy;  Matl_Gxz (x, 1 )  =  Gxz ; Matl_Gyz (x, 1 )  =  Gyz; 

%Second  ply 

[ ~ , Ex, Ey, ~ , ~ , Vxy, Vyx, Gxy, Gxz ,  Gyz ,  ~ ,  ~ ,  ~ ,  ~ ]  =  .  .  . 

LAMINATE_NASA_HALPIN_HYBRID_BASELINE_FUNC (40E- 
3, b (x) , h, alpha, theta (1,2)  ,  0 )  ; 

Matl  Ex(x,2)  =  Ex;  Matl  Ey(x,2)  =  Ey;  Matl  Vxy(x,2)  =  Vxy; 

Matl^Vyx (x, 2 )  =  Vyx; 

Matl_Gxy (x, 2 )  =  Gxy;  Matl_Gxz (x, 2 )  =  Gxz ; Matl_Gyz (x, 2 )  =  Gyz; 

%Third  ply 

[ ~ , Ex, Ey, ~ , ~ , Vxy, Vyx, Gxy, Gxz , Gyz  ,  ~ ,  ~ ,  ~ ,  ~ ]  =  .  .  . 

LAMINATE_NASA_HALPIN_HYBRID_BASELINE_FUNC (4  0E- 
3, b (x) , h, alpha, theta (1,3)  ,  0 )  ; 

Matl_Ex(x,3)  =  Ex;  Matl  Ey(x,3)  =  Ey;  Matl  Vxy(x,3)  =  Vxy; 

Matl^Vyx (x, 3 )  =  Vyx; 
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Matl_Gxy (x, 3 )  =  Gxy;  Matl_Gxz (x, 3 )  =  Gxz ; Matl_Gyz (x, 3 )  =  Gyz; 

%  Use  largest  value  of  Ex  to  decide  value  for  the  effective  modulus 
%  This  prevents  the  width  of  some  sections  from  being  too  large. 

Eff  Mod  =  max (Matl  Ex (x,  :));  %determines  max  Ex 

width (x,l)=  Matl_Ex (x, 1 ) /Ef f  Mod;  %makes  the  width  of  each  section 
a  scaler,  with  largest  being  =  1 

width (x, 2) =  Matl_Ex (x, 2 ) /Ef f  Mod; 
width (x, 3) =  Matl_Ex (x, 3 ) /Ef f  Mod; 

%  Prevent  I  (12) *2  from  being  >  I  (11) +1(22) 
if  width(x,l)<  width (x, 2 ) /2 . 5 
width (x, 1 ) =width (x, 2 ) ; 

h  new=h/ (3- (width (x, 2 ) / (min (width (x,  :))))); 
h=h  new; 

end 

if  width(x,3)<  width (x, 2 ) /2 . 5 
width (x, 3 ) =width (x, 2 ) ; 

h  new=hl / ( 3- (width (x, 2 ) / (min (width (x,  :))))); 
h=h  new; 

end 

%  Calculate  the  effective  density  for  the  Beam  Sections 
real  area  =  hl*b(x);  %real  cross  sectional  area  of  the  element 
faux_area  =  h*b (x) * (width (x, 1) /3  +  width (x, 2) /3  +  width (x, 3) /3) ; 
%area  of  I-beam  section  for  the  element 

area_scale  =  real_area/faux_area;  %scale  factor  to  multiply  density 
by  so  mass  is  same  for  all  elements 
Density (x) =den_cf *area_scale; 

%  prints  out  in  (1,  h,  bl (base  width) ,  b2  (top  width) ,  tl  (base) , 
t2 (top) ,  t3 (mid  width) ) 

fprintf (fid2, ' %  #6.12g,  %  #6.12g,  %  #6.12g,  %  #6.12g,  %  #6.12g,  % 

#6 . 12g,  %  #6 . 12g\n ' , . . . 

h/2 , h, b (x) * width (x, 3 ) , b (x) * width (x, 1 ) ,h/3,h/3,b(x) * width (x, 2 ) ) ; 
%prints  out  the  orientation  of  the  beam,  same  for  all  cases 
fprintf ( fid2 ,' 0 ., 1 ., 0 .  \n'); 

%  put  thickness  back  to  original  thickness 
h=hl  ; 

end 

%  This  section  is  for  the  Quad  (S8R)  elements 

fprintf (fid2,  ' **  Section:  Section-%-2 . Of \n ' ,  ( set_ctr+l ) ) ;  %  create 
section/profile 

fprintf (fid2, ' *Shell  Section,  elset=  PickedSet%-2 . Of , 
material=Material-%-2 . Of \n  '  .  .  . 

, (set_ctr+l) , (set_ctr+l) ) ; 
fprintf (fid2,  %6.12g,  %d\n ' , shell_thick, 5 ) ; 

%  This  section  is  for  the  TRI  (STRI65)  elements 

fprintf (fid2,  ' **  Section:  Section-%-2 . Of\n  ,  (set_ctr+2) ) ;  %  create 
section/profile 

fprintf (fid2, ' *Shell  Section,  elset=  PickedSet%-2 . Of , 
material=Material-%-2 . Of \n  '  .  .  . 


155 


, (set_ctr+2) , (set_ctr+2) ) ; 
fprintf (fid2, ' %6 . 12g,  %d\n ' , shell_thick, 5) ; 

%%  End  Instance 

fprintf ( f id2 , ' *End  Instance\n ' ) ; 
fprintf ( f id2 , ' **\n ' ) ; 

%%  Nset  for  Boundry  Conditions  and  Point  Load 

%Create  a  set  of  nodes  where  the  boundry  conditions  will  be  applied 
fprintf (fid2, ' *Nset,  nset=  PickedSet%-2 . Of ,  internal,  instance=Part- 1- 

1 \n '  .  .  . 

, (set_ctr+3) ) ; 

%Put  in  the  nodes  you  wish  to  clamp  in  this  case  the  base 
fprintf (fid2, ' %d, \n 1 , Base_Node_Loc) ; 

%Creat  a  set  for  the  Pt  Load  to  be  applied 

fprintf (fid2, ' *Nset,  nset=  PickedSet%-2 . Of ,  internal,  instance=Part- 1- 

1 \n '  .  .  . 

, ( set_ctr+4 ) ) ; 

%Put  in  the  node  where  you  want  the  force  applied, 
fprintf (fid2, ' %d, \n ' , Pt_Node_Loc) ; 

%%  End  Assembly 

fprintf ( fid2 , ' *End  AssemblyVn ' ) ;  %ends  the  assembly  portion  of  it 
fprintf ( f id2 , ' **\n ' ) ; 

%%  Materials 
%Standard  Heading 
fprintf (fid2, ' **  MATERIALS \n ' ) ; 
fprintf ( f id2 , ' **\n ' ) ; 

%  Determine  and  Assign  Matl  Prop  for  the  Veins  (B32  elements) 
for  x=l : size (angle,  1 ) 

theta  =  angle (x, 2 ) +theta_original ( 1 , : )  ; 

%set  up  material  properties 

fprintf (fid2, ' *Material,  name=Material-%-2 . Of \n ' , x) ; 

fprintf ( f id2 , ' *Density\n ' ) ; 

fprintf (fid2, ' %#6 . 6f\n ' , Density (x) ) ; 

fprintf (fid2,  ' *Elastic,  type=ENGINEERING  CONSTANTS \n ' ) ; 

%  Use  the  material  propertie  saved  from  Section  section 
%  determine  which  properties  to  use 
ctr=max (Matl_Ex (x,  : ) )  ; 
if  Matl_Ex(x,l)  ==  ctr; 
ctr  =  1; 

elseif  Matl_Ex(x,2)  ==  ctr; 
ctr=2 ; 

elseif  Matl_Ex(x,3)  ==  ctr; 
ctr=3 ; 

end 

fprintf (fid2, ' %#12 . 6g,  %#12.6g,  %#12.6g,  %#12.6g,  %#12.6g,  %#12.6g, 
%#12.6g,  %#12 . 6g, \n ' . . . 

,  Matl_Ex (x, ctr ) , Matl_Ey (x, ctr ) , Matl_Ey (x, ctr) ,0.1, 0.1, 0.1,  .  .  . 
Matl_Gxy (x, ctr) , Matl_Gxz (x, ctr) ) ; 
fprintf (fid2,  %#12.6g,\n', Matl_Gyz (x, ctr ) ) ; 
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end 


%  Assign  the  Matl  Prop  to  the  Quad  elements  (S8R) 

fprintf ( f id2  ,  ' *Material,  name=Material-%-2 . Of \n ' ,  (set  ctr+1)); 

fprintf ( f id2 , ' *Density\n ' ) ; 

fprintf ( f id2 , '%#6.6f\n', den_sh) ; 

fprintf (fid2, '*Elastic\n') ; 

fprintf (fid2,  ' %#12 . 6g,  %#6 . 6f \n ' ,  [shell_E  shell_v] ) ; 

%  Assign  the  Matl  Prop  to  the  TRI  elements  (STRI65) 

fprintf ( fid2 *Material,  name=Material-%-2 . Of \n ' , (set  ctr+2) ) ; 

fprintf ( f id2 , ' *Density\n ' ) ; 

fprintf (fid2, ' %#6 . 6f \n ' , den_sh) ; 

fprintf (fid2, ' *Elastic\n') ; 

fprintf ( f id2 ,  '%#12.6g,%#6.6f\n',  [ shell_E  shell_v] ) ; 

%%  Steps 

%create  the  steps  to  output  the  Eignevector/Eigenvalue  and  Point  load 
%response 

%%  Step  1  Point  Load 

%Create  a  static  step  to  evaluate  a  point  load  located  on  the  wing 
%Standard  Heading 

fprintf (fid2 , ' **  - 

- \n') ; 

fprintf ( f id2 , ' **  \n ' ) ; 

fprintf (fid2, ' **  STEP:  Step-l\n'); 

fprintf ( f id2 , ' **\n ' ) ; 

%Since  we  want  a  static  load  it  will  look  like  this 

fprintf (fid2, ' *Step,  name=Step-l,  perturbation\n ' ) ;  %  perturbation  must 
be  set  otherwise  step  will  carry  over  loads 
fprintf (fid2, ' *Static\n ' ) ; 


%Boundry  Conditions 
%Standard  Heading 
fprintf (fid2,  ' **  \n  '  )  ; 

fprintf (fid2, ' **  BOUNDARY  CONDITIONS\n ' ) ; 
fprintf (fid2,  ' **  \n ' )  ; 


%create  a  boundry  condition  (in  this  case  it  will  be  clamped) 
fprintf ( fid2 ,'* *  Name:  BC-1  Type:  Displacement/Rotation\n ' ) ;  %since 
only  one,  can  fix  6  dof,  says  type  of  boundry  cond 

fprintf ( fid2 *Boundary\n ');  %tells  abaqus  there  is  a  boundry  conditon 
fprintf ( fid2  ,  PickedSet%-2 . Of ,  1,  l\n',  (set_ctr+3) 
dir 

fprintf (fid2 , '^PickedSet%-2 . Of ,  2,  2\n', (set_ctr+3) 
dir 

fprintf (fid2, ' _PickedSet%-2 . Of ,  3,  3\n ' , ( set_ctr+3 ) 
dir 

fprintf (fid2 , '^PickedSet%-2 . Of ,  4,  4\n', (set_ctr+3) 
about  1  axis 

fprintf ( fid2 , 1 _PickedSet%-2 . Of ,  5,  5\n ' , ( set_ctr+3 ) 
about  2  axis 


;  %hold  in  x  or  1 

;  %hold  in  y  or  2 

;  %hold  in  z  or  3 

;  %hold  in  rot 

;  %hold  in  rot 


157 


fprintf (fid2, ' _PickedSet%-2 . Of ,  6,  6\n  , (set_ctr+3) ) ;  %hold  in  rot 
about  3  axis 

%  Loads 
%Heading 

fprintf ( f id2 , ' **\n ' ) ; 
fprintf (fid2, ' **  LOADS\n'); 
fprintf ( f id2 , ' **\n ' ) ; 

%Apply  they  type  of  load,  a  concentrated  node 

fprintf ( fid2  ,'* *  Name:  Load-1  Type:  Concentrated  f orce\n ' ) ; %name  of 
load 

fprintf ( fid2 *Cload\n '); %  type  =  concentrated 

fprintf (fid2, ' _PickedSet%-2 . Of ,  3,  %d\n ' , (set_ctr+4 ) ,  Pt_Load) ; 

%  Outputs 
%Standard  Heading 
fprintf ( f id2 , ' **\n ' ) ; 

fprintf (fid2, ' **  OUTPUT  REQUESTS \n ' ) ; 
fprintf ( f id2 , ' **\n ' ) ; 

%Displacement/Stress  outputs  (in  this  case  it  does  not  change) 
fprintf (fid2, ' **  FIELD  OUTPUT:  F-Output- 1 \n ' ) ; 
fprintf (f id2 , ' **  \n '  )  ; 

fprintf ( fid2 *Output,  field,  variable=PRESELECT\n ' ) ; 

%History  Output 
fprintf ( f id2  ,  ' **\n ' )  ; 

fprintf (fid2,  ' **  HISTORY  OUTPUT:  H-Output-1 \n ' ) ; 
fprintf ( f id2 , ' **\n ' ) ; 

fprintf (fid2, ' *Output,  history,  variable=PRESELECT\n ' ) ; 

%  End  Step 

fprintf ( fid2 ,' *End  Step\n'); 

%%  Step  2  Frequency 

%Create  a  step  that  wil  levaluate  Eigenvectors/Eigenvalues  of  the  wing 
%Standard  Heading 

fprintf ( fid2  ,  '**  - 

- \n') ; 

fprintf (fid2,  ' **  \n ' )  ; 

fprintf (fid2, ' **  STEP:  Step-2\n'); 

fprintf ( f id2 , ' **\n ' ) ; 

%Since  we  want  to  find  freqency,  it'll  look  like  this 
fprintf (fid2, ' *Step,  name=Step-2,  perturbation\n'); 
fprintf ( fid2 ,' Find  Freq.  \n'); 

fprintf ( fid2 ,' *Frequency,  eigensolver=Lanczos ,  acoustic  coupling=on, 

normalization=displacement\n ' ) ; 

fprintf (fid2,  ' %6  .  Of ,  ,  ,  ,  , \n ' , num_modes ) ; 

%Boundry  Conditions 
%Standard  Heading 
fprintf (fid2,  ' **  \n ' )  ; 

fprintf (fid2, ' **  BOUNDARY  CONDITIONS\n ' ) ; 
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fprintf (fid2, ' ** 


\n ' )  ; 


%create  a  boundry  condition  (in  this  case  it  will  be  clamped) 
fprintf ( fid2 ,'* *  Name:  BC-1  Type:  Displacement/Rotation\n ' ) ;  %since 
only  one,  can  fix  6  dof,  says  type  of  boundry  cond 

fprintf ( fid2 *Boundary\n ');  %tells  abaqus  there  is  a  boundry  conditon 
fprintf ( fid2  ,  1  PickedSet%-2 . Of ,  1,  l\n',  (set_ctr+3) 
dir 

fprintf ( fid2 , :  PickedSet%-2 . Of ,  2,  2\n', (set_ctr+3) 
dir 

fprintf (fid2, ' _PickedSet%-2 . Of ,  3,  3\n ' ,  ( set_ctr+3 ) 


;  %hold  in  x  or  1 
;  %hold  in  y  or  2 


dir 

fprintf ( fid2  ,  ' ^PickedSet%-2 . Of ,  4,  4\n',  (set_ctr+3) 
about  1  axis 

fprintf (fid2, ' _PickedSet%-2 . Of ,  5,  5\n ' , ( set_ctr+3 ) 
about  2  axis 

fprintf (fid2,  ' _PickedSet%-2 . Of ,  6,  6\n ' ,  (set_ctr+3) 

about  3  axis 


;  %hold  in  z  or  3 
;  %hold  in  rot 
;  %hold  in  rot 
;  %hold  in  rot 


%Outputs 

%Standard  Heading 
fprintf ( f id2  ,  ' **\n '  )  ; 

fprintf (fid2, ' **  OUTPUT  REQUESTS \n ' ) ; 
fprintf ( f id2 , ' **\n ' ) ; 


%Frequency  outputs  (in  this  case  it  does  not  change) 
fprintf (fid2, ' *Restart,  write,  f requency=0 \n ' ) ; 
fprintf (fid2, '  **  \n  '  )  ; 

fprintf (fid2, ' **  FIELD  OUTPUT:  F-Output-2 \n ' ) ; 
fprintf (fid2,  '  **  \n  '  )  ; 

fprintf ( fid2 ,' *Output ,  field,  variable=PRESELECT\n ' ) ; 

%End  First  Step 

fprintf ( fid2 ,' *End  Step\n'); 

%%  Close  the  File 

fclose ( f id2 ) ; 


%%  Run 

%  This  line  runs  the  completed  file 
%Make  matlab  slow  down  so  Abaqus  doesn't  crash 
pause  ( 5 ) 

eval ( [ ' dos ( ' ' abq6111  job='  Input_File  '  interactive'')']) 

%%  Open  . dat  File 

%Gather  the  results  from  the  .dat  file  and  store  in  a  matrix 
fid3=f open ([ Input  File ( 1 : length ( Input  File) -4)  ' . dat ' ] , ' r+ ' ) ;  %  Opens 

.dat  file  for  reading/writing 
idx  =  0 ; 

while  ~feof(fid3) 
idx  =  idx  +  1; 
tline=fgetl (fid3) ; 

if  isempty (strfind (tline, 'EIGENVALUE  OUTPUT  ')) 
==  0  %  Scan  till  *Nodes 
for  i=l:5 
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tline=f getl ( f id3 ) ;  %  Read  5  lines  of  blank  stuff 

end 

for  i=l : 1  %  for  i  =  1  to  the  number  of  modes  you  wish  to  get 
%Get  the  1st  Modal  Freq. 

Eigenvalue  string  ( i ,:)  =f getl  ( fid3 );  %Grabs  the 
eigenvector/value  date  from  . dat  file 

make_text  =  str2num  (Eigenvalue_string)  ;  %converts  text 
into  numbers  MATLAB  can  read 

Eigenvaluel  =  make  text  (1,4);  %  gets  the  Frequency  from 
the  string  of  numbers 

%Get  the  2nd  Modal  Freq. 

Eigenvalue  string  ( i ,:)  =f getl  ( fid3 );  %Grabs  the 
eigenvector/value  date  from  .dat  file 

make  text  =  str2num (Eigen  value_string) ;  %converts  text 
into  numbers  MATLAB  can  read 

Eigenvalue  =  make  text  (1,4);  %  gets  the  Frequency  from 
the  string  of  numbers 
end 

end 

end 

fclose ( f id3 ) ;  %  Closes  the  File 
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Appendix  H.  Wing  .inp  File 


Appendix  H  contains  a  condensed  version  of  the  .inp  file  submitted  to  ABAQUS 
in  order  to  solve  the  engineered  wing  FEA  Model.  Redundant  infonnation  has  been 
replaced  by  the  symbol  [snip] . . .  because  the  original  .inp  file  would  take  approximately 


658  pages  to  illustrate  due  to  the  vast  number  of  element  sets  and  material  properties. 


‘Heading 

“  Job  name:  BEAM_MODAL_BASELINE  Model  name:  08-Feb-2012 
*Preprint,  echo=NO,  model=NO,  history=NO,  contact=NO 

k  -k 

“  PARTS 

‘Part,  name=WING_FIT_SHOULDER_EE_08-Feb-2012 

■ k  • k 

‘End  Part 

■ k  • k 

““  ASSEMBLY 


‘Assembly,  name=Assembly 

k  : k 


‘Instance,  name=Part-l-l,  part=WING_FIT_SHOULDER_EE_08-Feb-2012 
*Node 


1,  0.0111069997, 

2,  0.0000000000, 

3,  0.0253340003, 

4,  0.0224740005, 

5,  0.0223999998, 


-0.0007809994,  0.0000000000 
0.0000000000,  0.0000000000 
0.0031570005,  0.0000000000 
0.0004820003,  0.0000000000 
-0.0012629998,  0.0000000000 


[snip] . 

22247, 

0.0405554337, 

0.0051652824, 

0.0000000000 

22248, 

0.0307590264, 

0.0053330138, 

0.0000000000 

k  k 

‘Element,  type=B32 
7356,  1105,  21460, 

15 

7357 , 

1106,  21477, 

1105 

[snip] . 

8447, 

380,  8626, 

381 

8448, 

381,  8579, 

7 

k  k 


‘Element,  type=S8R 


14, 

CO 
\ — 1 

19, 

1356, 

1363, 

7487, 

7488, 

7489,  7490 

15,  1356, 

1357, 

1502, 

1308, 

7491, 

7492, 

7493,  7494 

[snip] . 

7354, 

7433, 

7445, 

7446, 

7430, 

22002, 

21771, 

22234,  22229 

7355, 

7293, 

7244, 

7357, 

7292, 

22126, 

21749, 

22038,  21846 

k  k 
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*Element,  type=STRI65 

1,  1674,  1444,  1304,  7448,  7449,  7450 

2,  49,  2,  50,  7451,  7452,  7453 

[snip] ... 

6963,  7431,  7419,  7328,  21470,  21471,  21472 

6964,  7392,  7418,  7233,  21473,  21474,  21475 

-k  -k 

*Nset,  nset=  PickedSetl  ,  internal 
15,  1105,  21460 

*Elset,  elset=  PickedSetl  ,  internal 

7356 

*Nset,  nset=  PickedSet2  ,  internal 
1105,  1106,  21477 

*Elset,  elset=  PickedSet2  ,  internal 

7357 

[snip] ... 

*Nset,  nset=  PickedSetl093,  internal 
7,  381,  8579 

*Elset,  elset=  PickedSetl093,  internal 
8448 

*Nset,  nset=  PickedSetl094,  internal 


1,  3,  4, 

5, 

6,  7,  8, 

9, 

10, 

11, 

12,  13, 

14, 

15, 

16, 

17 

18,  19, 

20, 

21,  22, 

23, 

24, 

25, 

26,  27, 

28, 

29, 

30, 

31, 

32, 

33 

34,  35, 

36, 

37,  38, 

39, 

40, 

41, 

42,  43, 

44, 

45, 

46, 

47, 

48, 

49 

50,  51, 

52, 

53,  54, 

55, 

56, 

57, 

58,  59, 

60, 

61, 

62, 

63, 

64, 

65 

66,  67, 

68, 

69,  70, 

71, 

72, 

73, 

74,  75, 

76, 

77, 

78, 

79, 

80, 

81 

82,  83, 

84, 

85,  86, 

87, 

88, 

89, 

90,  91, 

92, 

93, 

94, 

95, 

96, 

97 

98,  99, 
112,  113 

100, 

101,  102, 

103, 

104, 

105,  106, 

107, 

108, 

109,  110, 

111, 

[snip] ... 

22238,  22239,  22240,  22241, 
22248, 

■ k  • k 

*Elset,  elset=  PickedSetl094 

22242,  22243,  22244, 

,  internal 

22245, 

22246, 

22247, 

14,  15, 

16, 

17,  18, 

19, 

20, 

21, 

22,  23, 

24, 

25, 

26, 

27, 

28, 

29 

30,  31, 

32, 

33,  34, 

35, 

36, 

37, 

38,  39, 

40, 

41, 

42, 

43, 

44, 

45 

46,  47, 

48, 

49,  50, 

51, 

52, 

53, 

54,  55, 

56, 

57, 

58, 

59, 

60, 

61 

62,  63, 

64, 

65,  66, 

67, 

68, 

69, 

70,  71, 

72, 

73, 

74, 

75, 

76, 

77 

78,  79, 

80, 

81,  82, 

83, 

84, 

85, 

86,  87, 

88, 

89, 

90, 

91, 

92, 

93 

94,  95, 
109 

96, 

97,  98, 

99, 

100, 

101 

,  102, 

103, 

104, 

105 

,  106, 

107, 

108, 

110,  111 

,  112,  113, 

114 

,  115 

,  116,  117, 

118 

,  119 

,  120, 

121, 

122 

,  123, 

124, 


[snip] ... 

7341,  7342,  7343,  7344,  7345,  7346,  7347,  7348,  7349,  7350,  7351,  7352, 
7353,  7354,  7355, 

■ k  ~k 

*Nset,  nset=  PickedSetl095,  internal 

1,  2,  3,  4,  5,  6,  7,  8,  9,  10,  11,  12,  13,  524,  525,  526 
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527,  528,  529,  530,  531, 
1153, 

532, 

533, 

534,  535,  536 

,  537,  538, 

539, 

540, 

[snip] ... 

6549,  6550,  6551 

6552,  6553,  6955,  6956, 

k  k 

6957, 

6958, 

6959,  6960, 

6961,  6962, 

6963, 

6964, 

*Elset,  elset=  PickedSetl095, 
1,  2,  3,  4,  5,  6,  7,  8,  9,  10, 
527,  528,  529,  530,  531,  532, 
1153, 

internal 

11,  12,  13,  524, 
533,  534,  535,  536 

525,  526 
,  537,  538, 

539, 

540, 

[snip] ... 

6552,  6553,  6955,  6956, 

6957, 

6958, 

6959,  6960, 

6961,  6962, 

6963, 

6964, 

k  k 


**  Section:  Section-1  Profile:  Profile-1 
‘Beam  Section,  elset=  PickedSetl  ,  material=Material- 1  , 
temperature=GRADIENTS ,  section=I 
3.75000000000e-005,  7 . 50000000000e-005,  0.000200000000000, 

0.000200000000000,  2.50000000000e-005,  2 . 50000000000e-005, 

0.000200000000000 

0.,1.,0. 

**  Section:  Section-2  Profile:  Profile-2 
[snip] ... 

**  Section:  Section-1093  Profile:  Profile-1093 
‘Beam  Section,  elset=  PickedSetl093,  material=Material-1093, 
temperature=GRADIENTS ,  section=I 
7.50000000000e-005,  0.000150000000000,  0.000150000000000, 

0.000150000000000,  5.00000000000e-005,  5 . 00000000000e-005, 

0.000125086683873 

o.,i.,o. 

**  Section:  Section-1094 

*Shell  Section,  elset=  PickedSetl094,  material=Material-1094 
2e-005,  5 

**  Section:  Section-1095 

*Shell  Section,  elset=  PickedSetl095,  material=Material-1095 
2e-005,  5 
*End  Instance 

k  k 

*Nset,  nset=  PickedSetl096,  internal,  instance=Part-l-l 

2, 

48, 

51, 

280, 

283, 

942, 

1269, 

*Nset,  nset=  PickedSetl097,  internal,  instance=Part-l-l 
1108, 

*End  Assembly 

k  k 

**  MATERIALS 
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*Material,  name=Material-l 

*Density 

3150.400000 

*Elastic,  type=ENGINEERING  CONSTANTS 

2 . 09090e+010,  6 . 92733e+009,  6 . 92733e+009,  0.100000, 

0.100000,  5 . 15227e+009,  5 . 15227e+009, 

3 . 78346e+009, 

*Material,  name=Material-2 
[snip] ... 

*Material,  name=Material-1093 

*Density 

1667.518713 

*Elastic,  type=ENGINEERING  CONSTANTS 

1 . 13487e+010,  9 . 51111e+009,  9 . 51111e+009,  0.100000, 

0.100000,  5 . 60787e+009,  5 . 60787e+009, 

5 . 09636e+009, 

*Material,  name=Material-1094 

*Density 

1249.600000 

*Elastic 

2 . 50000e+009,  0.300000 

*Material,  name=Material-1095 

*Density 

1249.600000 

*Elastic 

2. 50000e+009, 0.300000 

k  k  _ 


**  STEP:  Step-1 

k  k 

*Step,  name=Step-l,  perturbation 
*Static 


**  BOUNDARY  CONDITIONS 


**  Name:  BC-1  Type:  Displacement/Rotation 
*Boundary 

_PickedSetl096,  1,  1 
_PickedSetl096,  2,  2 
_PickedSetl096,  3,  3 
_PickedSetl096,  4,  4 
_PickedSetl096,  5,  5 
_PickedSetl096,  6,  6 

k  k 

**  LOADS 


**  Name:  Load-1  Type:  Concentrated  force 
*Cload 

_PickedSetl097,  3,  -1 . 000000e-003 

k  k 

**  OUTPUT  REQUESTS 

k  k 


**  FIELD  OUTPUT:  F-Output-1 

k  k 

*Output,  field,  variable=PRESELECT 


0.100000, 


0 . 100000, 
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k  k 


**  HISTORY  OUTPUT:  H-Output-1 

k  k 

*Output,  history,  variable=PRESELECT 
*End  Step 

k  k  _ 

k  k 

**  STEP:  Step-2 

k  k 

*Step,  name=Step-2,  perturbation 
Find  Freq. 

*Frequency,  eigensolver=Lanczos ,  acoustic  coupling=on, 
normalization=displacement 
9 

f  r  r  r  t 

k  k 

**  BOUNDARY  CONDITIONS 

k  k 

**  Name:  BC-1  Type:  Displacement/Rotation 
*Boundary 

_PickedSetl096,  1,  1 
_PickedSetl096,  2,  2 
_PickedSetl096,  3,  3 
_PickedSetl096,  4,  4 
_PickedSetl096,  5,  5 
_PickedSetl096,  6,  6 

k  k 

**  OUTPUT  REQUESTS 

k  k 

*Restart,  write,  frequency=0 

k  k 

**  FIELD  OUTPUT:  F-Output-2 

k  k 

*Output,  field,  variable=PRESELECT 
*End  Step 
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Appendix  I.  Wing  Multi-Run  Code 


This  code  was  used  to  create  multiple  runs  of  the  wing  code,  and  plot  the  results. 
This  Case  is  set  to  vary  the  laser  cut  angle. 

clc;  clear  all;  close  all; 

%%  Define  Number  of  Cases  to  Run 
NumCases=13; 

%%  Definite  Variation  in  Laser  to  Fiber  Orientation 
nstd  alpha=l; 
avg_alpha=0 ; 

alpha= [ -6; -5 ; -4 ; -3 ; -2 ; -1 ; 0 ; 1 ; 2 ; 3 ; 4 ; 5 ; 6] 
alpha= [alpha] ; 

%%  Define  Variation  in  Lamina  Layup  Orientation 
nstd_theta=l ; 
avg_theta=0 ; 

theta_original ( : , 2 ) =90-theta_original ( : , 2 ) ; 

%%  Define  Variation  in  Laminate  Thickness 
nstd_t=l ; 

avg_thickness=150 ; 
h  =  avg  thickness; 

h  =  h*lE-6;  %convert  from  integers  to  microns 
%%  Define  Parameter  Estimates 

[ahat, bhat,  ACI ,  BCI ]  =  unifit(theta  original ,.  0 1 )  ; 

%%  Run  the  Analysis 
for  t=l:NumCases 

[Eigen  valuel (t, : ) , Eigen  value2 (t, : ) ] =Final_Wing (alpha (t) , theta^origina 
1  ( t ,  :  )  ,  h  ( t )  ,  t )  ; 
end 

%%  Plot  the  Results  Histogram 
close  all 

%%  First  Modal  Frequency 
%  Open  a  Figure  the  size  of  the  screen 
scrsz  =  get (0, ' ScreenSize ' ) ; 

figure (' Position  1 ,  [1  scrsz (2)  scrsz (3)  scrsz (4) ] ) 

nbins=100 ; 
subplot  (331) 
hist (alpha, nbins) 
grid  on 

xlabel ([' Alpha  =  Mean  '  num2str(avg  alpha)  '  deg  |  std  =  ' 
num2str(nstd  alpha)]) 
ylabel ( ' Bins ' ) 

subplot  ( 332 ) 

hist (theta  original (:, 1 ), nbins ) 
grid  on 

xlabel ([' Theta  =  Mean  '  num2str(avg  theta)  '  deg  |  std  =  ' 
num2str (nstd_theta) ] ) 
ylabel ( ' Bins ' ) 
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subplot ( 333 ) 

hist (h ( : , 1) *1E6, nbins ) 

grid  on 

xlabel ([' Thickness  =  Mean  '  num2str(avg  thickness)  '  um  |  std  =  ' 
num2str  (nstd__t)  ]  ) 
ylabel ( ' Bins ' ) 

subplot ( 3 , 3 , 4 : 6 ) 

hist (Eigen  valuel , nbins ) 

grid  on 

xlabel ('1st  Modal  Frequency  (Hz)  ') 
ylabel ( ' Bins  '  ) 

title (['Mean  freq  =  '  num2str (mean (Eigen_valuel ) , ' %f ' )  '  (Hz)']) 

subplot ( 3 , 3 , 7 : 9 ) 

plot (1 :NumCases, Eigen_valuel , '  .  '  , 'MarkerSize ' , 8) 
grid  on 

xlabel ( ' Sample ' ) 
ylabel (' Frequency  (Hz) ') 
hold  on 

plot([0  NumCases ], mean (Eigen  valuel)*[l  1 ] , ' --g ; , ' LineWidth ' , 2 ) 
legend ( ' Samples  , ' Mean ' , ' Location ' , ' SouthEast ' ) 

title (['Mean  freq  =  '  num2str (mean (Eigen_valuel ) , ' %f ' )  '  (Hz),  Max  freq 

=  '  num2str (max (Eigen_valuel )  ,  ' %f ' )  '  (Hz),  Min  freq  =  ' 

num2str (min (Eigen^valuel )  ,  %f ' )  '  (Hz)']) 

hold  off 


%%  2nd  Modal  Frequency 
scrsz  =  get (0, ' ScreenSize ' ) ; 

figure ('Position', [1  scrsz (2)  scrsz (3)  scrsz (4) ] ) 

nbins=100 ; 
subplot  (331) 
hist (alpha, nbins) 
grid  on 

xlabel ([’ Alpha  =  Mean  '  num2str(avg  alpha)  '  deg  |  std  =  ' 
num2str(nstd  alpha)]) 
ylabel ( ' Bins ' ) 

subplot  ( 332 ) 

hist (theta  original (:, 1 ), nbins ) 
grid  on 

xlabel ([' Theta  =  Mean  '  num2str (avg^theta)  '  deg  |  std  =  ' 
num2str (nstd_theta) ]  ) 
ylabel ( ' Bins ' ) 

subplot ( 333 ) 
hist(h(:,l)*lE6, nbins) 
grid  on 

xlabel ([' Thickness  =  Mean  '  num2str(avg  thickness)  '  um  |  std  = 
num2str (nstd_t) ] ) 
ylabel ( ' Bins ' ) 

subplot ( 3 , 3 , 4 : 6 ) 
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hist  (Eigen  value2 , nbins ) 
grid  on 

xlabel('2nd  Modal  Frequency  (Hz)  ') 
ylabel ( ' Bins ' ) 

title (['Mean  freq  =  '  num2str (mean (Eigen_value2 ) ,  ' %f  '  )  '  (Hz)']) 

subplot ( 3 , 3 , 7 : 9 ) 

plot ( 1 : NumCases , Eigen  value2,  , ' MarkerSize ' , 8 ) 
grid  on 

xlabel ( ' Sample ' ) 
ylabel (' Frequency  (Hz)') 
hold  on 

plot([0  NumCases ], mean (Eigen  value2)*[l  l],'--g  ,  LineWidth ' , 2 ) 
legend ( ' Samples  1 ,  ' Mean '  ,  'Location  :  ,  ' SouthEast ' ) 

title (['Mean  freq  =  '  num2str (mean (Eigen_value2 )  ,  ' %f '  )  '  (Hz),  Max  freq 

=  '  num2str (max (Eigen^value2 )  ,  1 %f ' )  '  (Hz),  Min  freq  =  ' 

num2str (min (Eigen_value2 )  ,  ' %f ' )  '  (Hz)']) 

hold  off 
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Appendix  J.  Additional  Ideal  Wing  Mode  Shape  Images 


First  mode  shape: 
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Second  mode  shape 


a 
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Third  mode  shape 


V 
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Fourth  mode  shape: 
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